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Abstract 


The  Semi-spectral  Primit’ve  Equation  Model  (SPEM),  authored  by  Dr.  Dale  Haidvogel  of 
the  Johns  Hopkins  University,  is  one  approach  to  regional  and  basin  scale  ocean  modeling. 
This  user’s  manual  for  SPEM  describes  the  model  equations  as  well  as  what  the  user  must 
do  to  configure  the  model  for  a  specific  application. 
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Chapter  1 
Introduction 


This  user’s  manual  for  the  semi-spectrai  primitive  equation  model  (SPEM)  attempts  to 
describe  the  model  equations  as  well  as  what  the  user  has  to  do  to  configure  the  model  for  a 
specific  application.  Some  initial  tests  A  the  SPEM  are  described  in  Haidvogel,  Wilkin  and 
Young  [8]. 

The  principle  attributes  of  the  model  are  as  follows: 

General 

•  Primitive  equations  with  temperature,  salinity  and  an  equation  of  state. 

•  Hydrostatic  and  Boussinesq  approximations. 

•  Second-moment  conservation. 

•  Optional  Lagrangian  floats. 

Horizontal 

•  Orthogonal-curvilinear  coordinates. 

•  Arakawa  C  grid. 

•  Periodic  channel  or  user  supplied  boundary  conditions. 

Vertical 

•  Sigma  (bottom- topography  following)  coordinate. 

•  Chebyshev  modal  expansion. 

•  Rigid  lid  on  top. 

Chapters  2  ar'.'l  3  describe  the  model  physics  and  numerical  techniques  and  are  an  update 
of  the  description  in  Haidvogel,  Wilkin,  and  Young[S].  Chapter  4  lists  the  model  subroutines, 
functions  and  variables.  As  distributed.  SPEM  is  ready  to  run  a  seamount  resonance  prob¬ 
lem.  The  process  of  configuring  SPEM  for  a  particular  application  is  described  in  chapter 
5,  including  a  discussion  ot  the  seam  aunt  application.  Section  6  describes  the  Lagrangian 
floats  and  how  to  use  them,  including  the  plotting  post-processor  fltplt. 

The  SPEM  uses  version  3.0  of  the  NCAR  graphics  package  to  plot  the  model  fields  during 
execution.  You  must  have  version  3.0  if  you  want  to  use  the  SPEM  plotting  subroutines  as 
they  are.  You  will  have  to  edit  the  subroutine  velvet  as  described  in  section  4.5  in  order 
to  get  the  velocity  vectors  to  plot.  In  the  main  part  of  SPEM,  the  graphics  is  isolated  to 
calls  to  opngks.  peplt  (twice),  ploth  and  clsgks.  and  could  be  commented  out  if  desired. 


1.1  Acquiring  the  SPEM  code 

The  version  of  the  model  described  in  this  document  is  available  over  Internet  via  anonymous 
ftp  (user  ftp,  any  password)  from  jupiter.ino.ucar.edu  (IP  number  128.160.2.21).  Once  you 
are  connected,  type  ‘cd  spemS.O’  and  then  'dir  to  get  a  directory  listing.  You  will  find  the 
files: 


Readme 

Read  this  first 

Makefile 

For  a  UNIX  environment 

spem.com 

For  a  VMS  environment 

spem.  input 

The  input  text  hie 

spem3.0. whole 

The  model  code  without  include  directives 

spem3.0.f 

The  mode!  code  with  include  directives 

spem3.0.h 

The  SPEM  include  file 

mud2.f 

The  multigrid  solver 

mud2.doc 

The  documentation  for  mud2 

saxpy.f 

Cray  library  routines 

ezgrid.f 

A  program  for  making  a  simple  rectangular  grid 

fltplt. whole 

A  program  for  making  color  plots  from  the  float 
history  files  (without  includes) 

fltplt.f 

A  program  for  making  color  plots  from  the  float 
history  files  (with  includes) 

fltplt. h 

The  fltplt  include  file 

grid 

A  subdirectory  which  has  the  grid-generation 
programs  described  in  Hedstrom  and  Wilkin[10] 

mud2  is  copyrighted  by  NCAR  and  they  request  that  you  not  make  any  changes  to  it  without 
their  permission. 

If  you  are  unable  to  acquire  the  code  in  this  fashion  feel  free  to  contact  me  at: 

Kate  Hedstrom 

Institute  for  Naval  Oceanography 
Building  1103 

Stennis  Space  Center,  MS  39529 
(G0D-688-5737 

Internet:  kate^ncar. ucar.edu 
SPAN:  N  S  FG  VV : :  near .  ucar .  ed  u !  kate’’ 

1.2  Future  improvements 

A  number  of  improvements  to  the  model  are  being  worked  on  or  need  to  be  tested  more 
thoroughly  before  being  released.  These  include: 


•  A  first  and  second  moment  diagnostics  package. 


•  A  plotting  post-processor. 

•  A  netCDF  interface.  It  will  will  be  used  for  porting  binary  data  files  from  one  machine 
to  another  and  as  an  interface  to  the  ECMOP  facility  at  INO. 

•  A  surface  mixed-layer  in  SPEM.  One  option  which  has  been  investigated  is  the  Price- 
YVeller-Pinkel  mixed  layer  model. 

•  A  data  assimilation  scheme  in  SPEM. 

•  A  program  for  generating  polynomials  other  than  the  Chebyshevs.  Although  alternate 
polynomials  do  very  well  in  a  flat-bottomed  configuration,  we  are  still  exploring  their 
behavior  in  the  presence  of  topography. 

•  A  free  sea  surface. 

•  A  method  for  masking  out  land  areas. 

•  Multitasking  on  4-16  Cray  processors. 


1.3  Changes  between  version  2.0  and  version  3.0 

Some  quite  substantial  changes  have  been  made  to  the  model,  including: 

•  By  default,  the  SPEM  conserves  second  moments  to  within  a  few  percent  when  the 
energy-containing  scales  of  motion  are  well  resolved.  There  is  now  an  option  to  conserve 
second  moments  exactly. 

•  The  old  implementation  of  the  horizontal  pressure  gradient  terms  only  works  for  moder¬ 
ate  bottom  slopes  (2  x  10-2  maximum).  Some  corrections  can  be  added  to  the  pressure 
gradients  to  more  accurately  deal  with  steep  topography.  However,  these  corrections 
are  not  guaranteed  to  solve  all  problems  with  steep  topography. 

•  Separate  equations  for  temperature  and  salinity  and  an  equation  of  state  are  now 
used  instead  of  the  density  equation.  If  density  is  a  function  of  temperature  only,  the 
calculation  of  salinity  can  be  turned  off. 

•  There  is  an  option  for  calculating  the  two  components  of  the  surface  pressure  gradient. 

•  I'he  option  of  vertically  mixing  regions  of  static  instability  has  been  added. 

•  The  vertical  derivatives  and  integrals  using  model  polynomials  have  been  rewritten 
and  are  now  much  faster  on  vector  computers. 

•  Subroutines  for  following  Lagrangiun  floats  have  been  written.  There  is  also  a  program 
which  makes  color  plots  of  the  float  tracks. 

•  The  timestepping  has  been  changed  again.  See  Appendix  A. 


•  The  problem  which  SPEM  is  set  up  to  run  as  distributed  has  been  changed  from  a 
Kelvin  wave  to  a  seamount  resonance  problem. 

•  The  plotting  has  been  rewritten  to  take  advantage  of  the  new  features  of  NCAR  graph¬ 
ics  version  3.0,  especially  those  in  conpack. 

•  Several  changes  to  the  plotting  subroutines,  including  the  addition  of  slices  along  sur¬ 
faces  of  constant  rj. 

•  The  calls  to  mud2  have  been  modified  for  version  2.0  of  the  mud2  elliptic  solver. 

•  The  reading  and  writing  of  the  history  file  is  now  all  in  subroutine  histio. 

•  The  flag  slctpoly  has  been  taken  out  again  until  we  are  more  comfortable  with  using 
alternate  polynomials. 

•  The  common  blocks  have  been  put  in  an  include  file  and  alphabetized.  The  SPEM 
can  be  obtained  with  either  the  include  statements  and  the  include  file  or  with  the 
common  blocks  already  included. 

•  The  code  was  run  through  a  relabelling  program. 

1.4  Changes  between  version  1.0  and  version  2.0 

A  number  of  minor  changes  have  been  made  to  the  model,  including: 

•  The  timestepping  has  been  changed  from  leapfrog-trapezoidal  to  a  leapfrog  with  occa¬ 
sional  Adams-Bashforth  steps. 

•  A  flag  (perchan)  has  been  added  to  allow  the  domain  to  be  periodic  in  x. 

•  The  plotting  has  been  cleaned  up  and  the  plotting  subroutines  have  all  been  put 
together. 

•  The  model  density,  p,  has  been  redefined  so  that  ptotai  =  Po  +  p(~)  +  p{x,y,z,t). 

•  The  climatology  fields  uclm,  vclm,  and  rclm  have  been  added.  The  flags  uflags(9), 
vflags(9),  and  rflags(6)  switch  the  right-hand  side  terms  which  nudge  the  u,  v  and  p 
fields  towards  the  climatology  on  a  timescale  of  1  /rdmp. 

•  The  subroutine  getgrid  now  reads  in  h  and  f  as  well  as  the  grid  arrays.  The  program 
ezgrid  writes  out  a  rectangular  grid  with  its  topography  and  Coriolis  parameter. 

•  A  flag  (slctpoly)  has  been  added  to  allow  the  use  of  polynomials  other  than  the  Cheby- 
shevs.  These  polynomials  are  calculated  externally  by  a  program  which  is  available 
upon  request  and  which  will  be  documented  in  a  future  release. 

•  The  variables  fO,  beta,  xh,  yh,  dx.  dxsq,  de.  and  desq  have  all  been  deleted. 
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Chapter  2 

Model  Formulation 


2.1  Equations  of  motion 

The  primitive  equations  in  Cartesian  coordinates  can  be  written: 


du 

dt 

d<p  _  _ 

-f  v  ■  V  u  —  fv  =  -  —  +  Tu  +  Vu 

(2.1) 

dv 

dt 

+  V  ■  Vv  +  f  u  =  —  TT-  +  +  Vv 

dy 

(2.2) 

dT 

—  +  u  •  VT  =  Tt  +  Vt 
dt 

(2.3) 

P)Q 

™+t.vs  =  rs  +  Vs 

dt 

(2.4) 

P  =  p(T,S,P) 

(2.5) 

d<p  _  -pg 

dz  p0 

(2.6) 

du  dv  :  dw 
dx  ^  dy  ^  dz 

(2.7) 

where,  in  standard  notation: 

( u,v,w )  =  the  ( x,y,z )  components  of  vector  velocity  v 
p0  +  /?(x,  y,  z,  t)  =  total  density 
T{x,y,z,t)  =  total  temperature 
5(x,  ?/,  t)  =  total  salinity 
P  =  total  pressure  P  «  —p0gz 
4>(x,y,z,t)  =  dynamic  pressure  (pj pa) 
f(x,y)  =  Coriolis  parameter 
g  =  acceleration  of  gravity 

-  forcing  terms 


and 


(Vu,  V, ,,  Vt,  Vs)  =  diffusive  terms. 


Equations  (2.1)  and  (2.2)  express  the  momentum  balance  in  the  x  and  y  directions,  respec¬ 
tively.  The  time  evolution  of  the  perturbation  temperature  and  salinity  fields,  T(x,y,z,t) 
and  S(x,y,z,t),  are  governed  by  the  advective-diffusive  equations  (2.3)  and  (2.4).  The  equa¬ 
tion  of  state  is  given  by  equation  (2.5).  In  the  Boussinesq  approximation,  density  variations 
are  neglected  in  the  momentum  equations  except  in  their  contribution  to  the  buoyancy  force 
in  the  vertical  momentum  equation  (2.6).  Under  the  hydrostatic  approximation,  it  is  further 
assumed  that  the  vertical  pressure  gradient  balances  the  buoyancy  force.  Lastly,  equation 
(2.7)  expresses  the  continuity  equation  for  an  incompressible  fluid.  For  the  moment,  the 
effects  of  forcing  and  dissipation  will  be  represented  by  the  schematic  terms  T  and  V,  re¬ 
spectively.  The  horizontal  mixing  is  currently  implemented  as  either  harmonic  (Laplacian) 
friction,  biharmonic  friction  or  a  Shapiro  filter  (Shapiro  [13])  along  a  surfaces. 

2.2  Vertical  boundary  conditions 

The  vertical  boundary  conditions  can  be  prescribed  as  follows: 

top(z=0)  =  r£(x,y,t) 

vTz~  Ty(X^J^) 

«t| 7  =  T${x,y,t) 

7  =  Ts(x,y,t) 

w  =  0 

and  bottom  (z  =  —  h)  v  =  r^(x,  ?y,  t) 

uTz=  Ty  (X^V^) 

«r  §7  =  r£(x,y,t) 

«s§§  =  Ts(x,y,  t) 
v  ■  Vh  =  0. 

Surface  distributions  of  wind  stress  and  vertical  density  flux  (tj,t|)  are  prescribed 

on  the  top.  On  the  variable  bottom,  z  =  —h(x,y),  the  horizontal  velocity  components  are 
constrained  to  accommodate  a  prescribed  bottom  stress;  the  vertical  heat  and  salt  flux  is  also 
prescribed.  It  is  relevant  to  note  at  this  point  that  the  vertical  (polynomial)  discretization 
schemes  to  be  described  in  section  3.1  can  accommodate  quite  arbitrary  boundary  conditions 
on  a  =  0,  —h. 

2.3  Horizontal  boundary  conditions 

As  distributed,  the  model  configuration  is  that  of  a  periodic  channel.  If  the  variable  perchan 
is  set  to  false  then  the  boundary  conditions  are  that  of  a  closed  basin.  For  any  other  config¬ 
uration  it  is  the  responsibility  of  the  user  to  provide  the  appropriate  boundary  conditions  on 
u,  v ,  T,  5,  and  0,  where  is  the  transport  streamfunction.  At  every  time  step  the  subroutines 
bcs  and  bcpsi  are  called  to  fill  in  the  necessary  boundary  values. 
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If  the  biharmonic  friction  is  used,  a  higher  order  boundary  condition  must  also  be  pro¬ 
vided.  The  model  currently  has  this  built  into  the  code  where  the  biharmonic  terms  are 
calculated.  The  high  order  boundary  conditions  used  are  ^  (mn  §7%')  =  0  on  the  eastern 

and  western  boundaries  and  4-  (— )  =  0  on  the  northern  and  southern  boundaries.  The 

dy  \mn  dy 2  J 

boundary  conditions  for  u,T,  and  S  are  similar.  These  boundary  conditions  were  chosen 
because  they  preserve  the  property  of  no  gain  or  loss  of  momentum,  temperature,  or  salt 
across  a  closed  boundary. 


2.4  Sigma  (stretched  vertical)  coordinate  system 

From  the  point  of  view  of  the  computational  model,  it  is  highly  convenient  to  introduce  a 
stretched  vertical  coordinate  system  which  essentially  “flattens  out”  the  variable  bottom  at 
2  =  —h(x,  y).  Such  “sigma”  coordinate  systems  have  long  been  used,  with  slight  appropriate 
modification,  in  both  meteorology  and  oceanography  (e.g.,  Phillips  [11]  and  Freeman  et  al. 
[7]).  To  proceed,  we  make  the  coordinate  transformation: 


x  =  x 
V  =  V 

<7  =  1+  2  (z/h) 


and 


t  =  t. 


In  the  stretched  system,  the  vertical  coordinate  a  spans  the  range  —  1  <  <7  <  1;  we  are  there¬ 
fore  left  with  level  upper  (a  =  1)  and  lower  {a  =  —  1)  bounding  surfaces.  As  a  trade-off  for 
this  geometric  simplification,  the  dynamic  equations  become  somewhat  more  complicated. 


The  resulting  dynamic  equations  are,  dropping  the  carats 


and  remembering  the  chain  rule 


d± 

dx 


Z 


+  (1  -<r) 


hx  d <t> 
h"Wa 


du  .  _  _  d<f)  ,  ( gp  \  dh  ^  _ 

_-/0  +  ,.V»--s  +  (!-„)  + 

—  +  fu  +  vVv-~g-  +  (1~<7)  [—j-+F„  +  V, 

dT 

—  +  u-vr  =  Ft  +  vt 

at 

^  +  ?.VS  =  Ts  +  Vs 

at. 

P  =  p{T.  S,  P) 

do  /  - ghp\ 

da  \  2  p0  ) 

d  o  dCl 

—{hu)  +  -7-{lw)  +  k—  —  0 

dx  dy  da 


(2.8) 

(2.9) 

(2.10) 

(2.11) 

(2.12) 

(2.13) 

(2.14) 


where 


V  =  (u,u,fl) 

_  d  d  d 

v  •  V  =  u—  +  u—  +  ft— , 
ax  ay  aer 


and  the  vertical  velocity  in  sigma  coordinates  is 


=  -J- 


,1  d/i  .  d/i 

(1  -  <t)u  —  +  (1  -  a)V~Q^  +  2U; 


In  the  stretched  coordinate,  the  vertical  boundary  conditions  become: 

top(cr=  1)  (x)  fe  =  KiX’V't) 

I?  =  TT(x->y^) 

(nf)  ff  =  r|(x,y,<) 
ft  =  0 

(f)  If  =  rx  (*>?,*) 

(t)  ff  =  T>(x,y,t) 

(¥)U  =  o 

(¥)lf  =  0 

0  =  0. 


and  bottom  (<r  =  —1) 


Note  the  simplification  of  the  lower  (vertical  velocity)  boundary  condition  which  arises  from 
the  sigma  coordinate  transformation. 


2.5  Horizontal  curvilinear  coordinates 

In  many  applications  of  interest  (e.g.,  flow  adjacent  to  a  coastal  boundary),  the  fluid  may  be 
confined  horizontally  within  an  irregular  region.  In  such  problems,  a  horizontal  coordinate 
system  which  conforms  to  the  irregular  lateral  boundaries  is  advantageous.  It  is  often  also 
true  in  many  geophysical  problems  that  the  simulated  flow  fields  have  regions  of  enhanced 
structure  (e.g.,  boundary  currents  or  fronts)  which  occupy  a  relatively  small  fraction  of 
the  physical  /computational  domain.  In  these  problems,  added  efficiency  can  be  gained  by 
placing  more  (computational)  resolution  in  such  regions. 

The  requirement  for  a  boundary-following  coordinate  system  and  for  a  laterally  variable 
grid  resolution  can  both  be  met  (for  suitably  smooth  domains)  by  introducing  an  appropriate 
orthogonal  coordinate  transformation  in  the  horizontal.  Let  the  new  coordinates  be  £(x.y) 
and  y(x,  y )  where  the  relationship  of  horizontal  arc  length  to  the  differential  distance  is  given 
by: 

(ds)t  =  (— W  (2.15) 

V  »77  / 
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Ml  =  Q)  dr) 


(2.16) 


Here,  m(£,  tj)  and  n(£,7/)  are  the  scale  factors  which  relate  the  differential  distances  ( A£,  A p) 
to  the  actual  (physical)  arc  lengths. 

Denoting  the  velocity  components  in  the  new  coordinate  system  by 


v  ■  £  =  u 


(2.17) 


v  ■  fj  =  v 


(2.1S) 


the  equations  of  motion  (2.8)-(2.14)  can  be  re-written  (see,  e.g.,  Arakawa  and  Lamb  [4]): 

d  (  hu\  d  (  hu2\  d  (huv\  d  f  hutt\ 
dt  \mn  )  +  dp  \  m  )+  dcr  \  ran  J 


(2.19) 


d  (  hv\  d  f  huv\  d  ( hv2\  d  ( hvfi 

dt  J  y  n  mn 

+{{^)+vT((l)-urn  (£>}'m  = 

\m )  dr]  \2p0m /  dp 


(2.20) 


Ft  +  4  f  +  T  i~)  +  T  (~)  =:Ft  +  vT 

dt  \mn )  dc,  \  n  )  dp  \  m  J  dcr  \  mn  ) 


(2.21) 


4  (— )  +  4;  (— )  +  y  (— )  +  S-  (— )  =  jr,  +  p 

dt  \mn  J  d£  \  n  J  dp  \  m  J  dcr  \  mn  ) 


(2.22) 


p  =  p[T,S,P) 


(2.23) 


(2.24) 


o  ( ku\  d  ( hv\  d  I  hD. 
d(,  \  n  J  dp  \m  J  dcr  \mn 


All  boundary  conditions  remain  unchanged. 


(2.25) 
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Chapter  3 

Numerical  Solution  Technique 


3.1  Vertical  spectral  method 

The  vertical  (<r)  dependence  of  the  model  variables  is  represented  as  an  expansion  in  a  finite 
polynomial  basis  set  Pk((r)',  that  is,  we  set 

N 

btt,V,<T)  =  H  pk{<r)bk(^,v)  (3.1) 

t-=o 

where  the  arbitrary  variable  6,  introduced  here  for  convenience,  is  any  of  u,  v ,  p,  T,  S,  <p,  or  ft. 
The  only  restriction  placed  on  the  form  of  the  basis  polynomials  is  that  /ij  Pk{cr)da  —  £o,-, 
i.e.,  only  the  lowest  order  polynomial  may  have  a  non-zero  vertical  integral.  This  isolates 
the  external  mode  (or  depth  averaged  component  of  the  field)  in  the  amplitude  of  the  k  =  0 
polynomial.  In  practice,  the  spectral  technique  does  not  explicitly  solve  for  the  polynomial 
coefficients  b k  but  rather  for  the  actual  variable  values  at  “collocation'’  points  (or  equivalent 
grid  points)  crn  chosen  to  correspond  to  the  location  of  the  extrema  of  the  highest  order 
polynomial.  The  collocation  point  values  bn  are  thus  defined  as 


6(<rn)  =  £>,.(< Tn)bk 
k= o 


0  <  n  <  N 


and  will  be  functions  of  (£,77).  The  polynomial  coefficients  6^  can  be  recovered  from  the 
collocation  point  values  bn  by  a  linear  matrix  transformation.  Consider  a  matrix  F  whose 
elements  are 

Fnk  =  PkM  (3.3) 

and  let  b  and  b  be  column  vectors  whose  elements  are  bn  and  b respectively.  Then 


and  hence 


b  =  Fb 


b  =  F~xb. 


It  is  now  straightforward  to  represent  vertical  differentiation  and  integration  in  terms  of 
matrix  operators.  Consider  a  matrix  R  whose  elements  are 


Rnk  — 


dPk(<Tn) 


Then  the  matrix  Cpz  =  RF  1  will  perform  differentiation  of  the  model  field  (b)  in  the 
vertical  direction,  denoted  in  the  following  subsections  by  the  8a  operation,  since 

8ab  =  —k!fn)bk  =  Rb  =  RF~lb  =  CDZb.  (3.7) 

oa 


-L  =  Rb=  RF~xb  =  Cnzb. 


II 


Similarly,  consider  a  matrix  S  whose  elements  are 


Snk  =  f  Pk(cr)da-.  (3.8) 

Jan 

Then  the  matrix  C/jvr  =  SF~X  will  perform  vertical  integration,  denoted  in  the  following 
sections  by  the  J*  operation,  since 

rl  JL  . 

I]b=  I  bd(T  =  Y.  PM  da  bk  =  Sb  =  SF-'b  =  CiNTb.  (3.9) 

•*a  k= CT <7n 

By  default,  the  SPEM  uses  the  Chebyshev  polynomials.  The  matrix  operators  F.  R  and  S 
and  collocation  levels  an  for  this  basis  set  are  given  in  Appendix  B. 

3.2  The  collocation  method 

The  upper  and  lower  boundary  conditions  can  be  conveniently  implemented  for  the  poly¬ 
nomials  using  the  collocation  method.  The  boundary  conditions  are  of  the  form  where  the 
vertical  derivatives  of  the  fields  are  known,  for  instance 


This  comes  into  the  dynamical  equations  in  the  term 

d  (  du\ 

Tz  V  Tz  )  ' 

The  vertical  derivative  of  u  is  calculated  using  the  operator  Cqz  (equation  (3.7))  and  then 
the  values  of  i/|^  at  a  =  ±1  are  replaced  with  the  known  values  prior  to  the  second  vertical 
derivative  being  taken.  This  procedure  with  the  rules  for  vertical  integration  is  guaranteed 
to  preserve  the  global  balance  of  momentum,  heat,  and  salt. 

The  form  of  boundary  conditions  where  the  boundary  values  of  u  (or  v,  T  or  S)  are 
known  can  also  be  implemented  using  the  collocation  method.  In  that  case  the  dynamical 
equations  for  u  on  the  boundary  levels  are  replaced  by  setting  the  value  of  u  equal  to  the 
known  value.  This  approach  is  quite  commonly  used  in  finite-difference  techniques,  and  is 
extremely  convenient  when  the  boundary  conditions  provide  boundary  values  for  the  function 
being  sought. 

3.3  Horizontal  grid  system 

In  the  horizontal  (^,??),  a  traditional,  centered,  second-order  finite-difference  approximation 
is  adopted.  In  particular,  the  horizontal  arrangement  of  variables  is  as  shown  in  figure  3.1. 
This  is  equivalent  to  the  well-known  Arakawa  "C”  grid,  which  is  well  suited  for  problems 
with  horizontal  resolution  that  is  fine  compared  to  the  first  radius  of  deformation  (Arakawa 
and  Lamb  [4]), 


Figure  3.1:  Placement  of  variables  on  an  Arakawa  C  grid 

3.4  Conservation  properties 

SPEM  conserves  the  first  and  (optionally)  the  second  moments  of  u,v,S,  and  T,  where  the 
second  moments  of  velocity  correspond  to  kinetic  energy. 

3.4.1  First  moment  conservation 

To  conserve  first  moments,  the  four  vertical  advection  terms 

d_  d_  /MhA  d__  fhns\  and  d_  /hnr\ 

da  y  mn  )  da  \  mn  )  da  y  mn  J  da  \  mn  J 
are  computed  as  follows: 

d_  (httg\  _  fliQ. 
da  \  mn  )  \  mn 

where  q  denotes  any  of  the  four  quantities  u.  v.  S  and  T. 

To  ensure  the  conservation  of  the  first  moment  the  computed  advection  term  is  vertically 
integrated 

•AqBT  =  ^  J  Aq  da 

and  the  vertical  advection  term  is  corrected  by 

AqBC  —  Aq  —  AqBT- 


dq  (  hq  \  dQ 

-  ~T~  =  ^ 

da  \mn  J  da 


This  local  conservation  property  was  implemented  in  version  1.0  of  SPEM. 


3.4.2  Second  moment  conservation 

Version  3.0  will  also  conserve  second  moments  in  a  global  sense,  if  desired.  To  ensure  this 
conservation  of  a  quantity  <7,  the  following  integral  must  vanish  (this  can  be  seen  when  each 
term  is  integrated  by  parts  and  the  boundary  conditions  =  0  at  a  =  1,-1  are  used): 


Q°  =  £  (2 qA,BC  -  dc 


where  q  denotes  u,  v,  S  or  T .  The  horizontal  integral  of  these  deviations  from  zero 


Sq  =  JJ  Qd  d£  dr] 


is  projected  onto  the  existing  q  field  using 


Sq  =  JJJ  ^  (q  ■  qBc)  da  d£  dr) 


where  qBc  is  the  baroclinic  part  of  q  (qBc  =  q  —  |  /£/  qda).  The  correction 

.„  .  Sq 

AqBC  =  aiBC  +  -z-qBC 

leads  to  the  global  conservation  of  quadratic  moments  (i.e.  kinetic  energy,  T £  and  S'2). 

3.5  Semi-discrete  equations 

Using  the  representations  described  in  sections  3.1  and  3.3,  the  semi-discrete  form  of  the 
dynamic  equations  ( 2.19)— (2.25)  is: 


d  (  ul i"  \  c  c  (  uh 4 

Ot  I  m<n<  )  +  M  U  I  W 


+  6V  <un 


c 

+  M  — f—f 
min* 


m^v n 


-ci'nj-r.v  — 7-,  .  (.v  - 7  j  t  tv  1 

— — —  +  h*n6z  ( - )  -  u"  (  —  )  \v(-  = 

rn^n^  \nJ  \mJ 


—  +  ( 1  —  a )— — 

nK  2p0n< 


{  livTin  .  — 7TTiv  .  — TTT^l 

+  <  zzz—  +  T*  -  W  li^Sr,  (—)  1  TT>  = 

1  m<vn>n  *  \nJ  \mj  J 

T71  .  <jVip 

—  — -dy®  -4-  ( 1  —  a )— — — +  T>v  4-  T xi 
m'  z/?07?r' 


11 


hT\ 
mn  ) 

+  | 

'utfT*} 
,  "*  1 

>  +  | 

fvVT  1 
,  ^  J 

|  +  8  a  | 

f  hVtT 
[  mn 

'  hS^ 

)  +  «e 

( \ 

|  +  8V  < 

[+H 

f  has 

k  mn  y 

{  *  1 

l  ^  J 

[  mn 

=  Vf  +  J~T 


=  Vs  +  Ts 


p  -  plS,  T,  P) 


(3.12) 

(3.13) 

(3.14) 


6{a)  -  - 


gh 

,2  po 


l\p 


(3.15) 


+  4,{^}=0. 

mn 


(3.16) 


Here  6^  and  6n  denote  simple  centered  finite-difference  approximations  to  d/d £  and  d/dr). 

with  the  differences  taken  over  the  distances  A£  and  A 77,  respectively:  (  )  and  (  ) 

represent  averages  taken  over  the  distances  A£  and  At/;  8„  represents  a  vertical  deriva¬ 
tive  evaluated  according  to  the  prescription  described  in  Section  3.1;  and  1/  indicates  the 
analogous  vertical  integral.  The  vertical  advection  terms  are  actually  more  complicated,  as 
described  in  section  3.4. 

This  method  of  averaging  was  chosen  because  it  internally  conserves  mass  and  energy 
in  the  model  domain,  although  it  is  still  possible  to  exchange  mass  and  energy  through  the 
open  boundaries.  The  method  is  similar  to  that  used  in  \rakawa  and  Lamb  [4];  however, 
their  scheme  also  conserves  enstrcphy. 


3.6  Time-stepping:  depth-integrated  flow  (k  =  0  mode) 


In  continuous  form,  the  horizontal  momentum  and  continuity  equations  (2.19),  (2.20)  and 
(2.25)  can  be  written: 


d  (  hu\ 
dt  \  mn  ) 


Ru 


(  h  \  do 

\n)Ti 


(3.17) 


and 


d_  (hv\ 
dt  \  mn  ) 


d  ( hu\  d  I  hv\  d  (  hQ.  \  _ 

d(  \  n  )  dr)  \  m  )  da  \mn ) 


(3.18) 


(3.19) 


where  /?,,  and  /?„,  represent  the  sum  of  all  other  terms  in  the  u  and  v  equations,  respectively. 
Performing  a  vertical  average 


1 


l 


the  equations  become: 


^  —  =TTU-  ~ 


0  f  hn 
Ot  \mn 

0  /  hv 
8t  \mn 

d  ( hu 

dt\n 


h\  tty 
n)  dt 


(3.20) 


=  RV~  - 


h  \  tty 
m  J  dq 


d  ( hv\  ^ 
Or)  \m  ) 


where  the  overbar  indicates  a  vertically  averaged  quantity. 

By  virtue  of  (3.22),  the  depth- averaged  flow  is  horizontally  non-divergent.  This  enables 
us  to  represent  it  in  terms  of  a  streamfunction,  rb(t,q,t),  such  that 


U  =  -\~h 


n\  tty 
h)  dv 


(3.23) 


By  re-arrangement: 


f  m\  drb 

v  =  \TJ  ~dt 


1(1)  =  (T  Tmf 

dt  \m )  \h)  df 


4  (z) -(?)*-£• 

dt  \nJ  \  h  )  dr) 

Taking  the  curl  of  these  equations  yields  a  vorticity  equation  for  the  depth-averaged  flow: 
dq  did  fv\  d  fu\\  d  frn—\  d  (n—\_ 


(jls i  l  _  £.  (™tt\  -  —  (~ir\  -  r 

dt  \n)  dr)  \m)  J  d£\h  V  Or)  \h  V 


(3.25) 


or.  using  (3.23)  and  (3.24) 


m  \  d2v  /  n  \  d2rp  d  /  m  \  dv;  d  /  n  dv.' 
hn)  dt?  \hm)  dq2  +  dt,  \hn)  dt  +  dq  \  hrn )  dq 


The  introduction  of  the  streamfunction  i/>  serves  two  purposes.  First,  it  automatically  guar¬ 
antees  horizontal  non-divergence  of  the  transport  held,  as  required  by  (3.22).  Second,  it 
eliminates  any  dependence  on  the  depth -averaged  component  of  the  pressure  field  <t>  which 
contains  an  unknown  contribution  arising  due  to  the  rigid  lid. 

The  vorticity  equation  (3.26)  is  solved  by  first  obtaining  an  updated  value  of  q  by  appli¬ 
cation  of  the  leapfrog  (second-order)  time-differencing  scheme  described  in  Appendix  A.  The 
associated  value  of  the  streamfunction  is  determined  from  the  generalized  elliptic  equation: 


m  \  d2  V' 


hn)  dt? 


n  \  dzv  d  /  m  \  dv  0  /  n  \  ()v 

hrn)  dq2  +  dt  W in)  Ot  dq  \hrn  '  dq  ^ 


(3.27) 


A  solution  to  (3.27)  is  fully  prescribed  by  specifying  the  values  of  v  on  the  boundary  of  the 
model  domain.  The  SPEM  currently  uses  the  NCAR  elliptic  solver  mud2  \ see  Adams  [2]) 
to  solve  equation  (3.27). 


3.7  Time-stepping:  internal  velocity  modes  and  den¬ 
sity  field 

The  V  internal  modes  of  the  velocity  distribution  [(u,v)tJk  for  1  <  k  <  .V]  are  obtained 
by  direct  time-stepping  of  equations  (3.10)  and(3.11),  having  removed  their  depth- averaged 
component.  The  temperature  and  salinity  equations  (3.12)  and  (3.13)  can  be  similarly 
advanced  for  0  <  k  <  N.  See  Appendix  A  for  a  description  of  the  time-stepping  algorithms. 


3.8  Determination  of  the  vertical  velocity  and  density 
fields 

Having  obtained  a  complete  specification  of  the  u,  t>,  T,  and  S  fields  at  the  next  time  level 
by  the  methods  outlined  above,  the  vertical  velocity  and  density  fields  can  be  calculated. 
The  vertical  velocity  is  obtained  by  vertical  integration  (as  prescribed  in  section  3.1):  in 
particular: 

The  density  is  obtained  from  the  temperature  and  salinity  via  an  equation  of  state.  The 
SPEM  provides  a  choice  of  a  nonlinear  equation  of  state  p  —  p{T.S,z )  or  a  linear  equation 
of  state  p  —  p{T).  The  user  is  free  to  implement  their  preferred  equation  of  state. 

3.9  A  corrected  formulation  of  the  pressure  gradient 
terms 

The  pressure  gradient  terms  in  equations  (2.19)  and  (2.20)  are  written  in  the  form 

or 

hVo  +  ( 1  -  a)  (r~ , )  Vlu 

This  is  the  form  traditionally  used  in  sigma-coordinate  models  to  account  for  the  horizontal 
differences  being  taken  along  surfaces  of  constant  a.  This  form  can  be  shown  to  lead  to 
significant  errors  when  |V/i|  is  large  (Beckmann,  personal  communication). 

In  an  effort  to  reduce  these  errors  the  SPEM  will  optionally  calculate  additional  correction 
terms.  These  terms  are  based  on  a  Taylor  series  of  the  pressure  at  the  p  points  on  either  side 
of  a  velocity  point  (where  the  pressure  gradient  is  needed).  See  figure  3.2  for  the  geometry 
of  the  problem. 

The  pressure  values  which  are  required  are  pi(~n  —  Ac)  and  C>2(zn  +  Ac)  where  Ac  = 
l/l(<r  —  1  }  A  A .  This  can  be  expanded  in  c  coordinates  as: 

Oi  (cn  —  Ac)  —  Oi  (cn)  —  Ac  —  K 
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•essure  gradient  calculation 
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1  /  A  1 


,^l(2n) 


522 


_  1(A^)3 - 

6  ;  dz 3 


+ 


^(-n  +  A2)  =  4>2{Zn)  +  A2 


dfi^i  *-n  ) 


52 


,d2fi2{zn)  ,  1 


,d34>2(Zn) 


*  /  a  \2'-'  'V2\*-n)  r.  .3<-  r*v 

2(A2)  “a?-  +  6<a-’>  -&r-  + 


This  can  be  transformed  to  <7  coordinates  bv: 


/  2  \  dd>i(an) 

d>x(zn  -  A 2)  =  ^K)  -  A2  (j-J 

I,. 

'  UJ  ftr* 

02(^n  +  A2)  =  <f>2  (<7n)  +  A2 


2  \  ^^((Tn) 


2/  5(7 


2  \  2  d2fi2{crn) 


;) 


da2 


+  ■ 


The  horizontal  difference  of  these  terms  leads  to  the  following  form  of  the  pressure  gradient: 
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This  is  a  2-based  pressure  gradient  which  uses  the  vertical  differentiation  operator  (see  section 
•3.1)  to  extrapolate  to  the  depth  of  the  velocity  collocation  points.  If  both  o  points  are 
within  the  water  column  and  N  —  1  correction  terms  are  used,  this  calculation  is  equivalent 
to  summing  the  modes  to  find  6  at  the  given  depth.  Otherwise,  we  are  using  the  model 
polynomials  to  extrapolate  into  the  bottom. 


1!) 


3.10  Computation  of  the  surface  pressure 

In  timestepping  the  model  from  t  to  t  +  At,  most  of  the  terms  in  the  momentum  equations  are 
computed  directly.  However,  since  the  SPEM  has  a  rig'd  !id,  the  effects  of  the  surface  pressure 
gradient  are  not  immediately  known.  We  get  around  having  to  calculate  the  surface  pressure 
directly  by  calculating  the  curl  of  the  momentum  equations  and  finding  the  depth-integrated 
streamfunction  as  described  in  section  3.6.  However,  once  the  timestep  is  complete,  it  is 
possible  to  go  back  and  infer  what  the  surface  pressure  gradient  should  have  been. 


3.11  Convective  adjustment 

If  desired,  the  density  field  is  checked  for  occurences  of  static  instability  which  the  model 
will  then  attempt  to  remove.  It  cycles  through  the  water  column  three  times  using  a  simple 
mixing  scheme  which  conserves  T  and  S.  In  the  case  of  a  nonlinear  equation  of  state  or  of 
very  strong  advection  this  does  not  necessarily  guarantee  that  the  resulting  density  profile 
will  be  stable.  The  variables  T.  S  and  p  are  vertically  mixed,  but  the  corresponding  velocity 
fields  remain  unchanged. 

3.12  A  comment  on  timing  and  vectorization 

The  SPEM  is  written  in  highly  portable  Fortran  77  and  has  been  run  on  a  number  of  different 
computers.  A  version  with  42  x  41  gridpoints  and  7  polynomials  requires  the  amount  of  cpu 
time  per  timestep  shown  in  table  3.1.  This  version  has  no  diffusion  of  temperature  or  salinity 
and  does  not  conserve  second  moments;  turning  options  such  as  these  back  on  will  increase 
the  execution  time. 


Computer 

Time 

Precision 

Cray  Y-MP 

.083  sec 

64  bits 

Cray  X-MP 

.12  sec 

64  bits 

Cray  2 

.17  sec 

64  bits 

Stardent  Titan 

3.6  sec 

64  bits 

SGI  4D/220  (IRIS) 

2.3  sec 

32  bits 

Alliant  FX-80 

3.7  sec 

32  bits 

Alliant  FX-80 

6.4  sec 

64  bits 

Sun  SparcStation  1  + 

5.5  sec 

32  bits 

VAX  8800 

7.2  sec 

32  bits 

Table  3.1:  Relative  timings  for  several  different  computers.  Although  some  of  these  machines 
have  multiple  processors,  all  timings  were  done  on  one  processor 
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The  model  as  it  stands  has  been  optimized  for  the  Cray  and  uses  the  Cray  versions  of  the 
routines  saxpy  and  sdot  for  its  matrix  multiplication  and  dot  products.  If  one  were  to  run 
the  model  extensively  on  another  type  of  machine  it  may  be  worth  the  trouble  to  configure 
the  code  for  that  machine.  For  instance,  the  AUiant  has  a  compiler  directive  to  tell  it  that  a 
subroutine  call  within  a  loop  can  be  split  up  among  the  processors.  These  compiler  directives 
can  be  added  where  the  subroutine  call  is  independent  for  each  loop  index.  (However,  the 
most  obvious  of  these  loops  have  been  restructured  for  better  vectorization  and  are  no  longer 
independent.) 

On  supercomputers  such  as  the  Cray  X-MP,  efficient  operation  of  models  is  closely  tied  to 
the  degree  to  which  codes  can  be  vectorized.  In  the  present  SPEM,  it  is  important  to  realize 
that,  although  vertical  operations  such  as  integration  and  differentiation  axe  being  carried  out 
by  matrix  multiplication  (slow  transform)  methods,  such  operations  axe  nonetheless  highly 
vectorizable  when  being  simultaneously  applied  at  a  large  number  of  adjacent  horizontal 
points.  In  a  case  where  there  is  a  total  of  (say)  M  horizontal  gridpoints,  at  each  of  which  a 
vesical  operation  such  as  integration  is  needed,  the  requisite  M  integrations  can  be  obtained 
entirely  via  manipulations  of  vectors  of  length  M .  Since  in  most  applications  M.  will  oe  many 
hundreds  or  thousands,  and  will  typically  be  orders  of  magnitude  greater  than  the  number  of 
vertical  expansion  functions  (N),  vectorization  of  vertical  operations  by  horizontal  gridpoint 
is  extremely  powerful.  This  model  takes  full  advantage  of  this  vectorization  technique. 


Chapter  4 

Details  of  the  Code 


4.1  Main  subroutines 

A  flow  chart  for  the  main  program  is  shown  in  figure  4.1.  The  boxes  refer  to  subroutines 
which  are  described  as  follows: 

bcs  Horizontal  boundary  conditions  are  imposed  upon  the  u,v,T,S  and  ip  fields.  This 
can  be  modified  by  the  user  to  provide  solid  boundaries  or  open  boundaries  of  your 
favorite  type.  If  the  perchan  switch  is  true  then  the  boundaries  are  periodic  in  £. 

init  Does  everything  that  needs  to  be  done  to  start  up  the  model  run.  It  reads  initial 
parameters  and  u,  v,T,  S,  and  ip  fields  from  disk  or  calls  peminit.  It  then  calculates 
the  remaining  initial  fields  and  sends  out  the  initial  plots  and  opens  the  restart  file. 
The  flow  chart  for  init  is  shown  in  figure  4.2. 

omega  Calculates  the  vertical  velocity. 

prsgrd  Calculates  the  baroclinic  pressure  gradients,  with  an  entry  for  calculating  surface 
pressure  gradients. 

rhocal  Calculates  p  using  the  equation  of  state. 

srhs  Calculates,  stores,  and  tabulates  contributions  to  the  right  hand  side  of  equation 
(2.22),  where  the  advective  terms  have  been  moved  to  the  right  hand  side. 

tmstp  Performs  the  time  step  on  T,  S  and  the  time  step  on  the  ‘baroclinic  part'  of  u  and  v 
via  entry  points  tstp,  sstp.  ustp.  and  vstp  respectively.  The  entry  vrtstp  performs 
the  time  step  on  £,  solves  for  ip,  and  completes  the  time  step  on  u  and  v. 

trhs  Calculates,  stores,  and  tabulates  contributions  to  the  right  hand  side  of  equation 
(2.21),  where  the  advective  terms  have  been  moved  to  the  right  hand  side. 

urhs  Calculates,  stores,  and  tabulates  contributions  to  the  right  hand  side  of  equation 

(2.19) ,  where  all  the  terms  other  than  have  been  moved  to  the  right  hand  side. 

vrhs  Calculates,  stores,  and  tabulates  contributions  to  the  right  hand  side  of  equation 

(2.20) ,  where  all  the  terms  other  than  j^(  —  )  have  been  moved  to  the  right  hand  side. 

vrtrhs  Calculates  the  R.q  term  in  equation  (3.26). 
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4.2  Other  subroutines  and  functions 

Initialization 


chbset  Sets  up  Chebyshev  polynomial  functions. 

getgrid  Reads  in  the  curvilinear  coordinate  arrays  as  well  as  f  and  h  from  a  grid- 
generation  program. 

histio  Reads  or  writes  the  history  file  header. 

iput  Reads  model  fields  from  a  restart  record  on  disk, 
oput  Writes  model  fields  to  a  restart  record  on  disk. 

spoput  Writes  model  surface  pressure  gradients  to  a  restart  record  on  disk.  Sepa¬ 
rate  from  oput  so  that  the  initial  surface  pressure  gradient  can  be  calculated 
from  the  first  timestep. 

invmtx  Matrix  inversion  —  used  only  by  chbset. 
peminit  Sets  environmental  parameters  at  t  =  0. 

Numerical  techniques 

bndyc  Needed  for  mud2  but  never  used. 

coef  Specifies  the  coefficients  to  the  derivatives  in  mud2. 

dmean  Does  a  vectorized  vertical  integration  of  a  model  field  and  optionally  subtracts 
the  integral  (mean)  from  the  field. 

mud2  Generalized  elliptic  solver  (vectorized  for  Cray  computers;  Adams  [2]  and 
Adams  [3]).  Note:  mud2  only  needs  vorticity  values  in  the  interior  and  on  bound¬ 
aries  for  which  the  streamfunction  values  are  not  defined. 

rstmix  Mixes  vertically  in  regions  of  static  instability 

saxpy  Multiplies  a  constant  by  a  vector  and  adds  another  vector.  Used  in  multiplying 
matrices. 

scopy  Copies  a  vector  to  another  vector. 

sdot  Computes  an  inner  product  of  two  vectors. 

shapdl  Applies  a  one-dimensional  Shapiro  filter  to  a  vector. 

shapd2  Applies  a  two-dimensional  Shapiro  filter  to  a  matrix  by  calling  shapdl  in 
each  direction. 

sscal  Scales  a  vector  by  a  constant  scalar, 
ssum  Sums  the  elements  of  a  vector. 

Plotting  Warning:  for  the  slice  plots  an  array  which  is  dimensioned  L  x  M  x  0  :  N  is 
equivalenced  to  arrays  which  are  L  x  51  and  M  x  51.  This  will  lead  to  problems  for 
some  values  of  L.  M.  and  N.  Fix  this  by  commenting  out  the  equivalences  in  the 
subroutines  pltslca  and  pltslcb. 
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botplt  Plots  the  bottom  topography  in  slice  plots. 

cpmpxy  Replacement  for  the  conpack  mapping  routine  so  that  the  slice  and  slab 
contour  plots  come  out  in  the  appropriate  coordinate  system. 

cpshift  Makes  conpack  calls  in  such  a  way  as  to  emulate  our  old  modified  conrec.  It 
will  optionally  shift  the  contour  lines  by  half  a  contour  interval.  See  the  version 
3.0  NCAR  graphics  manual  [5], 

fill  Fills  an  array  for  plotting  horizontal  sections  of  the  model  fields  —  either  a 
level  of  constant  sigma  or  the  amplitude  of  a  vertical  mode. 

fx,  fy  Used  by  velvet  to  plot  the  horizontal  vectors  in  the  x,y  coordinates. 

getxxyy  Fills  the  xs,  ys,  xv,  yv,  xb,  and  yb  arrays  for  fx,  fy,  and  cpmpxy. 

label  Puts  labels  on  plots.  Changes  the  set  plotting  environment. 

mxf,  myf  Used  by  velvet  to  plot  vectors  in  x,y  coordinates. 

peplt  Creates  slab  plots  of  u,  v,  ip,  p ,  T,  S ,  w ,  £,  v,  <p  and  calls  pltslca  and  pltslcb  for 
the  slice  plots. 

ploth  Plots  bottom  topography,  its  gradient,  and  a  resolution  ratio  (see  Appendix 

C). 

pltperim  Draws  the  outline  of  the  model  domain  in  x,y  coordinates.  Used  with 

velvet. 

pltslca  Plot  slices  on  surfaces  of  constant  £. 

pltslcb  Plot  slices  on  surfaces  of  constant  //. 

velvet  A  modified  version  of  the  NCAR  Graphics  velocity  vector  subroutine,  allowing 
velocity  vectors  to  be  plotted  in  x,y  space  instead  of  £,77  space.  See  section  4.5 
and  Appendix  D. 

xyzplot  Fills  an  array  for  plotting  horizontal  sections  of  the  model  fields  on  a  level  of 
constant  z. 

yzplot  Fills  an  array  for  plotting  vertical  sections  of  the  model  fields  in  y,z  or  x,~ 
space. 

Lagrangian  floats 

bew  Boundary  conditions  on  w  (see  section  6.1.1). 

fltinit  Initializes  the  particle  positions,  called  from  the  main  program  after  init.  The 
particle  positions  should  be  specified  to  lie  in  the  range  1  <  £  <  L.  1  <  //  <  M. 
and  —1  <  a  <  1. 

fltstp  Called  by  the  main  program  to  advance  the  particles. 

rhoflt  Subroutine  which  returns  p  at  the  float  positions  (£,7/,<t)  by  using  the  equation 
of  state  p{T ,  S,  - ) . 

rk4  Fourth  order  Runge-Kutta  implementation  derived  from  Numerical  Recipes 
by  Press,  et  al.  [12].  Calls  xderiv,  yderiv.  and  zderiv  or  x2deriv,  y2deriv. 
and  z2deriv. 


sflt  Subroutine  which  returns  S'  at  the  float  positions  {£,r],cr)  by  a  bicubic  inter¬ 
polation  of  S  in  a  4  x  4  box. 

s2flt  Subroutine  which  returns  S  at  the  float  positions  77 ,  cr)  by  a  bilinear  inter¬ 

polation  of  S  in  a  2  x  2  box. 

tflt  Subroutine  which  returns  T  at  the  float  positions  (£,r],  a)  by  a  bicubic  inter¬ 
polation  of  T  in  a  4  x  4  box. 

t2flt  Subroutine  which  returns  T  at  the  float  positions  (£,7,  cr)  by  a  bilinear  inter¬ 
polation  of  T  in  a  2  x  2  box. 

xderiv  Subroutine  which  returns  d^jdt  at  the  float  positions  (<^,  77 .  cr )  by  a  bicubic 
interpolation  of  u  ■  m  in  a  4  x  4  box. 

x2deriv  Subroutine  which  returns  d£/dt  at  the  positions  (£,77,77)  by  a  bilinear  inter¬ 
polation  of  u  •  m  in  a  2  x  2  box. 

yderiv  Subroutine  which  returns  drj/dt  at  the  positions  (£,77,77)  by  a  bicubic  interpo¬ 
lation  of  v  •  n  in  a  4  x  4  box. 

y2deriv  Subroutine  which  returns  drj/dt  at  the  positions  (<J ,  77,  cr )  by  a  bilinear  inter¬ 
polation  of  v  •  n  in  a  2  x  2  box. 

zderiv  Subroutine  which  returns  dcr/dt  at  the  float  positions  (£,77,77)  by  a  bicubic 
interpolation  of  £1  in  a  4  x  4  box. 

z2deriv  Subroutine  which  returns  dcrjdt  at  the  positions  (£,77,77)  by  a  bilinear  inter¬ 
polation  of  0  in  a  2  x  2  box. 

Other 

vmax  Returns  the  largest  element  in  a  vector. 

vmin  Returns  the  smallest  element  in  a  vector. 

4.3  Important  parameters 

Following  is  a  list  of  the  important  parameters  in  the  model.  Most  other  parameters  are 

derived  from  L.  M,  and  N. 

L  Number  of  grid  points  in  the  £  direction. 

M  Number  of  grid  points  in  the  r\  direction. 

N  Number  of  baroclinic  polynomials  in  the  vertical. 

nflt  Number  of  Lagrangian  floats. 

nwrk  Size  of  the  work  space  used  by  mud2. 
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4.4  Common  blocks  and  the  variables  within  them 

All  variables  which  have  units  are  given  in  MKS  units. 

BLANK  The  main  time-dependent  model  fields 

time  time 

u  velocity  component  in  direction 

v  velocity  component  in  77  direction 

w  Q,  velocity  component  in  a  direction 

rho  perturbation  density  (total  density  =  pQ  -f  p) 
t  temperature  T 

s  salinity  S 

psi  horizontal  transport  streamfunction  ib 

vrt  vorticity  of  depth-averaged  flow  field  £ 

dwds  vertical  derivative  of  fi,  reallv  4-{— )  or  -MP- 

/chb/  Polynomial  transformation  matrices 

pi  7r  =  3.14159265  . . . 

sig  array  of  sigma  values  cr{k) 

cp  matrix  of  orthogonal  polynomials  F  =  Pk{crn)  (equation  3.3).  Multiply  a  mode 
amplitude  vector  by  cp  to  get  the  vector  of  values  at  the  an  locations.  b{crn). 

cf  inverse  of  the  matrix  of  orthogonal  polynomials.  F_1  (equation  3.5).  Multiply 

a  b(crn)  vector  by  cf  to  get  the  mode  amplitude  vector. 

cd  derivatives  of  Pj(cn)  (equation  3.6) 

cdz  derivative  transform  matrix  Cdz  (equation  3.7).  Multiply  a  vertical  vector  by 
cdz  to  get  the  vector  of  derivatives. 

cint  integral  transform  matrix  C/at  (equation  3.9).  Multiply  a  vertical  vector  bv 
cint  to  get  the  vector  o:  integr  j,  integrating  down  from  the  surface. 


csum 

integrals  of  Pk(cr )  over  —  1  <  a  <  +1 

/climat/ 

Model  climatology  used  in  nudging  terms 

tclm 

T  climatology 

sclm 

5  climatology 

uclm 

u  climatology 

vclm 

v  climatology 

/coord/  Orthogonal  coordiua'e  transformation  arrays 
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pm  coordinate  transformation  metric  m 
pa  am;  itinat e  transformation  metric  n 
dnd\  4-(~) 
dn.  ie  ■£-(—) 

at I'm  ' 

/ellip/  Parameters  needed  for  elliptic  soh-er.  See  the  mud2  documentation  for  details 

iprm  parameters  for  the  mud2  elliptic  solver 
ipnu  parameters  for  the  mud2  elliptic  solver 
ewrk  w.  rk  array  of  size  nwrk 

nmua  number  of  steps  between  printing  out  the  mud2  error 
mgopt  p  up  rioters  for  the  mud2  elliptic  solver 

/flags/  Logical  arrays  tor  the  light  hand  side  calculations 

uflags  logical  array  used  to  determine  which  terms  are  to  be  included  in  the  u  equa¬ 
tion 

vflags  logical  array  used  to  determine  which  terms  are  to  be  included  in  the  v  equa¬ 
tion 

sflags  logical  array  used  to  determine  which  terms  are  to  be  included  in  the  S  equa¬ 
tion 

tflags  logical  array  used  to  determine  which  terms  are  to  be  included  in  the  T  equa¬ 
tion 

mixflag  logical  flag  to  turn  on  vertical  mixing 

stpflag  logical  flag  to  determine  the  type  of  time  step  corrections  to  take,  true  is  for 
third-order  Adams- Bashforth  and  false  is  for  leapfrog-trapezoidal  steps. 

eqstflag  logical  flag  to  determine  which  equation  of  state  to  use.  true  for  a  non-linear 
p(T,S,z)  and  false  for  a  linear  p(T) 

/floats/  Lagrangian  float  variables 

fit  A  two-dimensional  array  containing  the  particle  information.  The  first  index 

specifies  which  particle  is  being  referred  to  and  has  a  range  of  1  to  nflt.  The 
second  index  specifies  the  field  according  to  the  key: 

1  =£ 

2  =  7 

3  =  <7 

4  =  =  um 

5  =  17 (£,»/, )  =  v  -n 

6  = 

7  -- 

S  =  S((,q,<r) 

!)  =  dUc-'/.<x' 
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kptime  Interval  between  disk  writes  of  float  information 

fltflag  Logical  flag  to  specify  whether  or  not  floats  are  being  used 

flt4flag  Logical  flag  to  specify  4  x  4  or  2  x  2  horizontal  interpolation,  true  for  4x4. 

nfltrrec  Number  of  float  restart  records  to  read  in  from  disk  during  initialization 

/grdpts/  x,y  (physical  space)  gridpoint  locations,  used  in  plotting  the  model  fields 

xp  x  coordinates  of  the  ip  points 

yp  y  coordinates  of  the  ip  points 

xmin  minimum  x  value  on  the  ip  grid 

ymin  minimum  y  value  on  the  ip  grid 

xmax  maximum  x  value  on  the  ip  grid 

ymax  maximum  y  value  on  the  ip  grid 

xl  length  of  domain  in  £  direction,  xl  =  xmax  —  xmin 

el  length  of  domain  in  y  direction,  el  =  ymax  —  ymin 

/id/  Descriptive  identification  for  the  current  simulation 

ident  160  character  string,  first  24  characters  are  used  to  label  plots 

/ionum/  Unit  numbers  for  I/O.  The  SPEM  does  not  specifically  open  files  or  use  any 
particular  file  names.  Fortran  compilers  have  default  file  names  for  the  unit  numbers, 
such  as  fort. 3  or  FOR001.DAT  and  these  are  the  file  names  used  by  SPEM  for  its 
binary  files. 

iin  unit  number  for  reading  in  history  file 

iout  unit  number  for  writing  out  history  file 

/parm/  Physical  constants  and  fields 
h  bottom  depth 

f  Coriolis  parameter  =  2  •  Q garth sin(0) 

uvnu2  lateral  Laplacian  mixing  coefficient  (constant)  for  u.v 
tnu2  lateral  Laplacian  mixing  coefficient  (constant)  for  T 
snu2  lateral  Laplacian  mixing  coefficient  (constant)  for  5 
uvnu4  lateral  biharmonic  mixing  coefficient  (constant)  for  u,  u 
tnu4  lateral  biharmonic  mixing  coefficient  (constant)  for  T 
snu4  lateral  bi harmonic  mixing  coefficient  (constant)  for  N 
rhoO  mean  density  pn 

g  acceleration  due  to  gravity 

rdrg  bottom  drag  coefficient 

:f  1 


hmin  minimum  depth  of  the  topography 
hmax  maximum  depth  of  the  topography 

dhdx 

dhde 

m  07) 

rfac  -3— 

2(>o 

rdmp  inverse  timescale  of  the  nudging  towards  climatology 
/perbcs /  Variables  particular  to  periodic  boundary  conditions 

perchan  logical  flag — true  for  periodic  boundary  conditions  in  the  x  direction 

perintg  logical  flag  used  only  in  the  periodic  channel  version — false  if  you  will  specify 
tp  on  both  channel  walls,  true  if  you  want  the  model  to  calculate  the  channel 
transport 

gfpy  used  to  calculate  along-channel  momentum  balance 

py  used  to  calculate  along-channel  momentum  balance 

grnfnc  the  transport  Green’s  function  for  a  periodic  channel 

/ pg/  Pressure  gradient  arrays  and  parameters 

phix  £  component  of  the  baroclinic  pressure  gradient 

phie  t)  component  of  the  baroclinic  pressure  gradient 

sphix  £  component  of  the  surface  pressure  gradient 

sphie  r]  component  of  the  surface  pressure  gradient 

npgc  number  of  pressure  gradient  correction  terms  to  calculate 

nspr  number  of  steps  between  calculating  the  surface  pressure  gradients 

/pit/  Plotting  flags  and  parameters 

nplvl  number  of  levels  to  be  plotted 

npl  array  of  specific  a  levels  to  be  plotted,  from  0  to  N 

nslcea  number  of  y-z  cross  sections  to  be  plotted 

nsla  array  of  specific  ^-locations  of  the  y-z  cross  sections 

nslceb  number  of  x-z  cross  sections  to  be  plotted 

nslb  array  of  specific  ^-locations  of  the  x-z  cross  sections 

plot  u  logical  switch  to  turn  on  plotting  of  the  u  field 

plot  v  logical  switch  to  turn  on  plotting  of  the  v  field 

plot  s  logical  switch  to  turn  on  plotting  of  the  5  field 

plot  t  logical  switcli  to  turn  on  plotting  of  the  T  field 

plot  w  logical  switch  to  turn  on  plotting  of  the  w  field 


plot  rho  logical  switch  to  turn  on  plotting  of  the  p  field 

plot  psi  logical  switch  to  turn  on  plotting  of  the  ip  field 

plot  vort  logical  switch  to  turn  on  plotting  of  the  £  field 

plot  vec  logical  switch  to  turn  on  plotting  of  the  v  field 

plot  spr  logical  switch  to  turn  on  plotting  of  the  surface  pressure  gradient 

zslab  determines  whether  the  level  (slab)  plots  are  of  planes  of  constant  a  or  constant 
a .  true  for  plots  at  a  constant  a 

zplot  array  of  specific  z  levels  to  be  plotted 

ztop  the  2- level  of  the  top  of  the  slice  plots 

zbot  the  2-level  of  the  bottom  of  the  slice  plots 

rbarsub  logical  flag  to  specify  whether  rhobar  is  subtracted  out  before  p  is  plotted 

/stp/  Time  step  parameters  (see  Appendix  A) 

dt  time  step  A t 

ntmes  number  of  steps  in  current  run 

nrst  number  of  steps  between  storage  of  restart  fields 

nrrec  number  of  restart  records  previously  written  to  disk 

nplt  number  of  steps  between  plots  (calls  to  subroutine  peplt) 

nmix  number  of  steps  between  vertical  mixings 

ncorr  number  of  steps  between  correction  steps 

iic  time  step  counter 

jjc  time  step  component  counter 

lnew  one  of  several  time  level  indices  (equal  to  either  1  or  2)  used  for  the  last  array 
index  in  the  model  fields,  lnew  points  to  the  current  time  level 

lold  pointer  to  the  previous  time  level 

Istp  pointer  to  the  time  level  to  which  the  current  changes  are  added 

ldis  pointer  to  the  time  level  used  in  the  dissipative  terms  (not  lrhs  for  numerical 
stability  reasons) 

lrhs  pointer  to  the  time  level  used  to  calculate  the  right-hand  side  terms 

Itrp  pointer  to  the  time  level  of  the  right  hand  side  arrays  (ru.  rv.  etc.)  to  which 

the  changes  are  being  added 

labl  pointer  to  a  time  level  for  the  right  hand  side  arrays 

lab2  pointer  to  a  time  level  for  the  right  hand  side  arrays 

cfactO  weighting  factor  for  the  right  hand  side  arrays 
cfactl  weighting  factor  for  the  right  hand  side  arrays 


cnb  weighting  factor  for  !F{t  —  A£) 

cnc  weighting  factor  for  T*‘*.  +  At) 

dts  size  of  the  timestep  between  lstp  and  lnew 

/user/  Arrays,  variables,  and  constants  defined  by  the  user  as  being  relevant  to  a  specific 
application 

/vec2/  velvet  internal  common  block  —  allows  the  plotting  of  a  velocity  vector  every  nth 
grid  point 

incx  plots  every  incx  vector  in  the  x  direction 
incy  plots  every  incy  vector  in  the  y  direction 

/vrtadv/  Needed  for  energy  conservation  operations 

nrgeons  logical  flag — true  if  you  want  to  conserve  energy 
c3d  3-d  work  array 

/wrk/  Temporary  work/storage  arrays 

rs  stores  right  hand  side  of  5  equation 

rt  stores  right  hand  side  of  T  equation 

ru  stores  right  hand  side  of  u  equation 

rv  stores  right  hand  side  of  v  equation 

rz  stores  right  hand  side  of  vorticity  equation 

tmp  2-d  work  array,  tmp  and  tmp2  are  used  to  carry  the  depth-averaged  forcing 
from  ustp  and  vstp  to  vorteq  and  should  not  be  changed  in  that  part  of  the 
main  loop 

tmp2  2-d  work  array 

wrz  1-d  (vertical)  work  array 

wvz  1-d  (vertical)  work  array 

a3d  3-d  work  array 

b3d  3-d  work  array 

/xxyys/  Temporary  x,  y  arrays  for  the  plots 

xs  x  position  for  fx  in  the  vertical  y ,  z  slice  plots 

ys  y  position  for  fy  in  the  vertical  y,  z  slice  plots 

xv  x  position  for  fx  in  the  horizontal  slab  plots 

yv  ij  position  for  fy  in  the  horizontal  slab  plots 

xb  x  position  for  fx  in  the  vertical  x.z  slice  plots 

yb  y  position  for  fy  in  the  vertical  x,  z  slice  plots 


Other  variables 


akuv  vertical  diffusive  coefficient  for  u,v  (statement  function  of 

akt  vertical  diffusive  coefficient  for  T  (statement  function  of  i,j,k) 

aks  vertical  diffusive  coefficient  for  5  (statement  function  of  i,j,k ) 

lspace  determines  whether  the  slab  plots  are  of  constant  cr  levels  or  the  mode  ampli¬ 
tude.  lspace  =  true  will  produce  plots  in  physical  space. 


4.5  Changes  to  velvet 

A  number  of  changes  have  to  be  made  to  the  NCAR  graphics  [6]  routine  velvet  in  order  for 
the  vector  plots  to  work  in  the  SPEM.  Appendix  D  shows  the  UNIX  'context  differences' 
between  the  unmodified  and  the  modified  versions  of  velvet. 

The  most  important  change  allows  the  plots  to  be  made  in  the  curvilinear  coordinates 
instead  of  rectangular  coordinates.  This  is  done  by  commenting  out  the  lines 

FX(X,Y)  =  X 
FY(X,Y)  =  Y 

and  the  lines 

MXF (XX , YY , UU , VV , SFXX , SFYY , MXX , MYY )  =  MXX  +  IFIX(SFXX*UU) 

MYF (XX , YY , UU , VV , SFXX , S FYY ,  MXX , MYY )  -  MYY  +  IFIX(SFYY*VV) 

and  providing  the  functions  fx,  fy,  mxf,  and  myf. 

Another  change  to  velvet  allows  the  vectors  to  go  outside  the  model  boundary.  Add  the 
lines 

call  gqclip(IER, ICLP , IAR) 
call  gsclip(O) 

before  the  line 

C  DRAW  THE  VECTORS 

This  will  save  the  current  clipping  state  and  disable  the  GKS  clipping.  The  two  lines 

CALL  GQCLIP(IER, ICLP, IAR) 

CALL  GSCLIP(O) 


a  few  lines  down  should  be  commented  out  so  that  the  saved  clipping  state  is  not  overwritten. 


Chapter  5 


Configuring  SPEM  for  a  Specific 
Application 


This  chapter  describes  the  parts  of  SPEM  for  which  the  user  is  responsible  when  configuring 
it  for  a  given  application.  Section  5.1  describes  the  process  in  a  generic  fashion  while  section 
5.2  steps  through  the  application  of  SPEM  to  a  seamount  resonance  problem.  As  distributed, 
SPEM  is  ready  to  run  the  seamount  example  in  a  periodic  channel.  Section  5.3  describes 
the  changes  which  are  needed  for  a  non-periodic  application. 


5.1  Configuring  SPEM 

5.1.1  Model  domain 

One  of  the  first  things  the  user  must  decide  is  how  many  grid  points  to  use  (and  can  be 
afforded).  There  are  three  parameters  which  specify  the  grid  size: 

L  Number  of  finite-difference  points  in  £. 

M  Number  of  finite-difference  points  in  tj. 

N  Number  of  baroclinic  vertical  modes. 

Since  the  model  uses  the  mud2  elliptic  solver  you  are  not  free  to  pick  just  any  values  for  L 
and  M.  For  a  non-periodic  domain  the  following  relations  must  be  satisfied: 

L  =  nxl  x  2(nsteP-1)  -f 1  (5.1) 

M  =  nyl  x  2(nsteP"1)  +  1  (5.2) 

where  the  solver  is  more  efficient  for  higher  values  of  nstep,  since  nstep  defines  the  number 
of  sub-grid  levels  used  in  the  multi-grid  algorithm.  If  periodic  boundary  conditions  are  being 
used  then  the  L  relation  becomes: 

L  =  nxl  x  2(nsteP~1)  +  2.  (5.3) 

See  the  mud2  documentation  for  a  description  of  the  mud2  solver. 

L,  M,  and  N  are  parameters  which  are  used  in  almost  every  subroutine.  They  should 

be  set  to  your  chosen  values  throughout  the  code.  This  is  also  a  good  time  to  initialize  the 

appropriate  parameters  for  mud2.  In  peminit  set 

iprm(6)  =  nxl 
iprm(7)  =  nyl 
iprm(8)  =  nstep. 


The  size  of  the  work  array  ewrk  is  given  by  the  parameter  nwrk.  For  periodic  domains 
nwrk  should  be  set  to  16  •  4  •  L  •  M/3  while  for  open  domains  the  size  14-4-LM/3is 
sufficient. 

5.1.2  x,y  grid 

The  subroutine  getgrid  is  called  by  peminit  to  read  in  the  grid  arrays,  the  topography, 
and  the  Coriolis  parameter.  The  file  it  reads  is  produced  by  the  grid  generation  programs 
described  in  Hedstrom  and  Wilkin  [lOj.  The  variables  which  are  read  by  the  subroutine 
getgrid  are: 

pm,  pn,  dndx,  dmde,  xp,  yp,  f,  h. 

From  these  fields  getgrid  calculates  the  related  variables: 

xmin,  xmax,  ymin,  ymax,  xl,  el,  hmin,  hmax,  dhdx,  dhde. 

The  aspect  ratio  of  the  domain  is  given  by  el/xl.  The  plotting  subroutines  currently  assume 
that  the  aspect  ratio  will  be  less  than  2.5.  If  this  is  not  the  case  then  you  can  adjust  the 
variables  xl,  x2,  yl  and  y2  in  the  subroutines  peplt  and  ploth  for  your  domain. 

For  a  simple  rectangular  geometry  a  short  program  called  ezgrid  can  be  used  to  calculate 
the  required  grid  arrays,  ezgrid  is  used  in  the  seamount  example  described  in  section  5.2. 

5.1.3  £,  7]  grid 

Before  providing  initial  conditions  and  boundary  conditions  the  user  must  understand  the 
model  grid.  The  fields  are  laid  out  on  an  Arakawa  C  grid  as  in  figure  3.1.  The  overall  grid 
is  shown  in  figure  5.1.  The  thick  outer  line  shows  the  position  of  the  model  boundary.  The 
points  inside  this  boundary  are  those  which  are  advanced  in  time  using  the  model  physics. 
The  points  on  the  boundary  and  those  on  the  outside  must  be  supplied  by  the  boundary 
conditions. 

The  three-dimensional  model  fields  are  carried  in  four- dimensional  arrays  where  the 
fourth  array  index  refers  to  the  two  time  levels  (see  appendix  A).  The  integers  i,  j,  and 
k  are  used  throughout  the  model  to  index  the  three  spatial  dimensions: 

i  Index  variable  for  the  f  direction, 

j  Index  variable  for  the  tj  direction. 

k  Index  variable  for  the  a  direction,  k  =  0  refers  to  the  bottom 

while  k  =  N  refers  to  the  surface. 

The  range  of  f  is  1  to  L  and  the  range  of  77  is  1  to  M.  Therefore  i  and  £  are  the  same  at 
points,  as  are  j  and  r]. 

5.1.4  Initial  conditions 

The  initial  values  for  the  model  fields  are  provided  by  peminit.  Only  the  interior  points  of 
u,u,T,  S  and  U’  need  be  initialized  in  peminit  because  bcs  will  be  called  after  peminit. 
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x  -  u  points 
□  -  v  points 
O  -  p  points 
•  -  xji  points 

Figure  5.1:  The  whole  grid 
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5.1.5  Equation  of  state 

The  equation  of  state  is  used  in  the  subroutines  rhocal  and  rhoflt.  Two  versions  are 
provided  in  SPEM  as  it  is  distributed:  a  nonlinear  p  =  p(T,S,z)  and  a  linear  p  =  p(T).  It 
is  possible  to  provide  your  own  equation  of  state,  but  if  you  use  floats  you  should  make  sure 
that  rhocal  and  rhoflt  are  consistent.  Also,  it  is  a  good  idea  to  write  it  in  a  vectorizable 
form  if  you  run  SPEM  on  a  vector  machine. 


5.1.6  rhobar  and  rhobarz 

The  statement  function  rhobar  is  used  in  prsgrd  to  designate  a  density  stratification  which 
is  a  function  of  2  only  and  which  will  be  subtracted  out  before  the  calculation  of  the  pressure 
gradients,  rhobar  is  also  set  in  some  of  the  plotting  routines  to  allow  the  plotting  of  p  —  ~p 
rather  than  the  plotting  of  p  alone.  There  is  a  related  statement  function  rhobarz  where  p 
is  expressed  as  a  function  of  2  as  opposed  to  i,  j,  and  k. 


5.1.7  Climatology 

The  climatology  fields  uclm,  vclm,  tclm  and  sclm  contain  the  values  to  which  u,v,T  and 
S  are  nudged  when  the  body  forcing  is  turned  on.  These  arrays  are  initialized  in  peminit. 
The  inverse  nudging  timescale  rdmp  is  also  initialized  in  peminit. 

5.1.8  Boundary  conditions 

The  horizontal  boundary  conditions  are  provided  by  the  subroutine  bcs  and  its  entry  bcpsi. 
They  are  called  for  every  timestep,  including  during  the  initialization,  and  must  provide  the 
boundary  values  for  the  fields  u,  v ,  T ,  S  and  ip.  The  user  configures  this  subroutine  to  provide 
the  boundary  conditions  desired. 

In  the  special  case  of  the  periodic  version  (perchan  =  true)  the  subroutine  bcs  sets  the 
boundary  values  to  wrap  around  appropriately.  Values  at  i  =  L  are  set  to  those  at  i  =  2, 
values  at  i  =  0  are  set  to  those  at  i  =  Lm2,  and  values  at  i  =  1  are  set  to  those  at  i  =  Lm. 
There  is  an  extra  set  of  p  points  and  u  points  for  compatibility  with  the  open  version. 

The  elliptic  solver  also  has  to  know  what  type  of  boundary  conditions  to  provide.  These 
are  specified  by  iprm(2)  through  iprm(5)  and  are  initialized  in  peminit. 


5.1.9  Timestep 

The  variable  dt  is  initialized  in  peminit.  For  simple  wave-like  problems  the  overall  stability 
is  limited  by  a  CFL-type  condition  such  that  c-^  <  1.  However,  when  horizontal  and  vertical 
diffusion  are  included,  the  stability  issue  becomes  more  complicated. 

It  is  sometimes  useful  to  be  able  to  change  the  timestep  and/or  the  horizontal  friction 
parameters  in  the  middle  of  a  model  run.  When  the  model  is  restarted  from  a  history  file, 
values  for  dt,  uvnu4,  etc.  are  read  in.  These  old  values  will  have  to  be  overwritten  in  init 
after  the  call  to  histio. 
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5.1.10  Model  Input 

ident  A  string  160  characters  long  used  to  identify  the  experiment  which  is  read  in  on 
unit  5.  The  first  24  characters  are  used  to  label  the  plots. 

5.1.11  Block  data 

As  many  of  the  user  choices  as  possible  are  put  into  the  block  data  at  the  end  of  the  program. 
They  are  as  follows  and  are  described  more  fully  in  section  4.4. 

sflags  Logical  flags  to  tell  the  model  which  terms  in  equation  2.22  to  use.  See  the  note  in 
the  block  data  tor  details. 

tflags  Logical  flags  to  tell  the  model  which  terms  in  equation  2.21  to  use. 

uflags  Logical  flags  to  tell  the  model  which  terms  in  equation  2.19  to  use. 

vflags  Logical  flags  to  tell  the  model  which  terms  in  equation  2.20  to  use. 

ntmes  Number  of  time  steps  to  take. 

stpflag  Logical  flag  to  specify  the  type  of  the  correction  timestep. 

ncorr  How  often  to  take  a  correction  step. 

nrrec  Number  of  restart  records  to  read  in  from  disk. 

mixflag  Logical  flag  which  is  true  if  vertical  mixing  of  statically  unstable  fields  is  desired, 
nmix  How  often  to  mix  vertically. 

eqstfiag  Logical  flag  which  is  true  if  the  non-linear  equation  of  state  is  desired, 
nplt,  nrst,  nmud  How  often  to  produce  various  types  of  output, 
g,  rhoO  Environmental  parameters. 

rdrg,  uvnu2,  snu2,  tnu2,  uvnu4,  snu4,  tnu4  Frictional  parameters. 

perchan,  perintg  Logical  flags  for  the  periodic  channel  boundary  conditions. 

nspr  How  often  to  calculate  surface  pressure  gradients. 

npgc  Number  of  corrections  to  the  pressure  gradient  terms  to  compute. 

nrgcons  Logical  flag  for  quadratic  moment  conservation. 

fltflag,  flt4flag,  nfltrrec,  kptime  Lagrangian  float  information. 

nplvl,  zslab,  npl,  zplot  Control  of  horizontal  level  plots. 

nslcea,  nsla,  nslceb,  nslb,  zbot,  ztop  Control  of  vertical  slice  plots. 

rbarsub  Logical  flag  which  is  true  if  p  —  ~p  should  be  plotted  instead  of  p. 

plot  . . .  Control  which  model  fields  are  plotted. 
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5.1.12  Vertical  diffusion 

The  vertical  diffusive  coefficients  akt,  aks,  and  akuv  can  vary  over  space.  They  are  currently 
implemented  as  statement  functions  within  trhs,  srhs,  urhs,  and  vrhs.  If  these  coefficients 
are  sufficiently  complicated  and  if  enough  memory  is  available,  it  is  sometimes  more  efficient 
to  implement  them  as  three-dimensional  arrays,  as  described  in  section  5.1.13. 

5.1.13  User  variables 

It  is  possible  for  the  user  to  add  new  variables  and  common  blocks  appropriate  to  a  given 
application,  for  instance,  some  boundary  conditions  require  that  you  keep  track  of  more 
information  than  the  model  currently  stores.  These  new  variables  can  be  initialized  in 
peminit  or  init,  but  recall  that  peminit  is  not  called  when  restarting  from  a  history  file. 
Variables  which  are  initialized  in  peminit  should  be  written  out  to  the  history  file,  either 
in  histio  (if  they  do  not  change),  or  in  oput. 

5.2  Seamount  Resonance  Example 

The  particular  application  for  which  we  have  configured  the  model  includes  a  diurnal  tide 
over  a  tall  seamount  in  a  stratified  ocean.  The  question  is  whether  or  not  the  tide  will  excite 
a  resonance  with  seamount-trapped  waves  which  have  a  period  close  to  one  day.  If  so,  this 
could  explain  the  large  bottom-trapped  currents  observed  over  Fieberling  Guyot. 

We  have  chosen  to  place  the  sea  mount  in  a  periodic  channel  to  simplify  the  boundary 
conditions.  The  grid  is  rectangular  with  the  finest  resolution  over  the  seamount  and  coarser 
resolution  elsewhere.  The  tidal  flow  is  provided  by  the  boundary  conditions,  specifically  the 
value  of  streamfunction  at  the  far  wall. 

5.2.1  Model  domain 

A  relatively  small  number  of  grid  points  were  chosen  for  reasons  of  economy.  Values  for 
L,  M,  and  N  are: 

L  =  42 
M  =  41 
N  =  7. 

These  values  lead  to  the  mud2  parameters  of 

iprm(6)  =  nxl  =  5 
iprm(7)  -  nyl  =  5 
iprm(S)  =  nstep  =  4. 

5.2.2  ezgrid 

For  this  geometry  one  has  a  choice  ot  using  the  grid-generation  programs  described  in  Hed- 
strom  and  Wilkin  [10]  or  of  creating  a  stand-alone  program  to  create  the  grid  file  analytically. 


A  short  program  called  ezgrid  was  written  to  produce  a  uniform  rectangular  grid  file.  This 
program  was  modified  to  produce  a  stretched  grid  with  a  Gaussian  seamount  when  the 
variable  stretch  is  true.  The  fluid  depth  is  5000  meters  and  the  seamount  is  4500  meters 
tall. 

To  use  ezgrid  for  another  application,  read  the  comments  in  the  code  and  make  the 
appropriate  changes.  For  a  grid  with  constant  grid  spacing,  set  stretch  to  false  and  supply 
the  size  of  the  domain  in  the  parameters  xl  and  el.  You  will  also  need  to  provide  it  with 
h(i,j)  and  f(i,j). 

5.2.3  Initial  conditions  and  the  equation  of  state 

We  would  like  the  initial  conditions  to  be  motionless  fluid  with  an  exponential  stratification. 
The  statement  functions  ufld,  vfld,  and  psifld  are  all  set  to  zero. 

The  stratification  can  be  provided  by  T ,  S,  or  both.  For  simplicity  we  will  only  have 
an  active  temperature  field  and  we  will  use  the  linear  equation  of  state  where  density  is  a 
function  of  temperature  only.  We  want  the  density  to  be  28.0  on  the  bottom  and  27.5  at  the 
top  with  an  e-folding  scale  of  1000  meters.  The  initial  temperature  is  set  to  1.43  +  3. 57e-/1000 
in  tfld. 

Since  density  does  not  depend  on  salinity,  we  have  a  choice  on  how  to  handle  the  salinity 
field.  We  can  either  use  it  as  a  passive  tracer  or  not  timestep  on  it  at  all  by  setting  sflags(7) 
to  false.  We  will  use  it  as  a  passive  tracer  and  initialize  it  in  sfld  to  be  a  function  of  y. 

5.2.4  rhobar  and  rhobarz 

We  have  set  rhobar  and  rhobarz  to  the  density  field  derived  from  the  initial  temperature 
field. 

5.2.5  Climatology 

We  have  turned  nudging  off  for  this  calculation  so  the  climatology  fields  are  not  used. 

5.2.6  Boundary  conditions 

The  boundary  conditions  are  walls  to  the  north  and  south  (rj  —  l.M).  The  walls  can  be 
implemented  simply  by  specifying  no  flow  out  of  the  box: 

u  =  0  &  j  —  l.M. 

In  the  presence  of  horizontal  diffusion  it  is  also  necessary  to  provide  horizontal  momentum 
and  density  fluxes  through  the  walls.  These  fluxes  are  determined  by  the  values  of  u  and 
p  outside  the  boundary.  No  flux  boundaries  were  chosen  by  setting  n0utside  =  ^inside  and 

^outside  =  Pinside- 

The  boundaries  to  the  east  and  west  are  periodic,  bcs  is  configured  to  provide  periodic 
boundary  conditions  when  perchan  is  true. 

Finally,  the  boundary  conditions  provide  the  barotropic  tidal  flow.  The  total  transport 
through  the  domain  is  given  by  the  difference  between  the  values  of  ip  on  the  northern  and 


southern  walls.  We  can  tell  SPEM  that  we  will  be  providing  the  value  of  ip  on  the  northern 
wall  by  setting  perintg  to  false  (the  southern  wall  is  always  set  to  zero). 

We  want  the  velocity  far  from  the  seamount  to  be 

u  —  i'0sm(^t) 

where  U0  is  1  mm/sec  and  ui  is  27r/ldav.  The  corresponding  value  of  can  be  gotten  by 
multip’ving  UQ  by  the  cross-sectional  area  of  the  channel  (5000  m  by  320  km).  This  leads 
to  $0  =  1.6  x  10b  or  1.6  Sverdrups.  We  do  not  want  to  shock  the  system  by  turning  on  the 
forcing  all  at  once  so  we  ramp  it  up  with  a 

5  [‘  + tanh  (^r). 

over  a  period  of  several  days. 

5.2.7  Timestep 

The  timestep  dt  is  initialized  in  peminit  to  be  a  day  over  12S. 

5.2.8  Block  data 

The  terms  to  be  used  in  the  equations  of  motion  are  specified  by  the  flags  in  the  block  data. 
We  want  the  Coriolis,  pressure  gradient,  advection  and  biharmonic  friction  in  the  momentum 
equations  so  uflags  and  vflags  2.  5,  7  and  8  are  true  while  the  rest  are  false.  The  value  of 
the  biharmonic  friction  coefficient  uvnu4  is  set  to  1  x  lO10.  We  want  an  adiabatic  calculation 
so  only  the  temperature  advection  tflags(5)  is  true.  To  enable  salinity  timestepping  and 
advection  we  set  sflags  5  and  7  to  true. 

The  variable  perchan  is  true  for  periodic  boundary  conditions  in  if,  but  set  perintg  to 
true  because  the  bcs  is  supplying  the  channel  transport, 
eqstflag  is  false  to  get  a  linear  equation  of  state. 

The  model  has  been  set  up  to  run  for  ten  days  or  1280  timesteps. 


5.2.9  Output 

The  model  writes  some  information  to  standard  oat: 

read  grid  file 

Stretched  grid  for  the  seamount  example 

Seamount  Example  Run  on  Cray  2 

May  15,  1990 

L  =42 

M  =41 

N  =  7 

xl  =  3 . 31 1E+05 

el  =  3 . 201E+05 
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rhoO  =  1.000E+03 
g  =  9 . 810E+00 
rdrg  =  3 . 000E-04 
hmin  =  5.000E+02 
hmax  =  5.000E+03 
dt  =  6 . 750E+02 
uvnu2=  1.000E+03 
snu2  =  1.000E+03 
tnu2  =  1.000E+03 
uvnu4=  1.000E+10 
snu4  =  1.000E+10 
tnu4  =  1.000E+10 
rdmp  =  0.000E+00 

initial  error  for  greens  function  solution  =  8.618E-05 
init  call  to  mud2:  ierr  =  -1 


Plotting  for  day:  0.00  has  been  done 


64 

3 . 426E-03 

128 

3 . 072E-03 

192 

2.719E-03 

256 

2.402E-03 

320 

2 . 289E-03 

320 

3 . 962E-06 

384 

2 . 203E-03 

448 

2 . 131E-03 

512 

2 . 068E-03 

576 

2 . 004E-03 

640 

1 . 936E-03 

640 

3 . 344E-06 

704 

1 . 862E-03 

768 

1.779E-03 

832 

1.694E-03 

896 

1 . 598E-03 

960 

1 .505E-03 

960 

2 . 591E-06 

1024 

1 . 410E-03 

1088 

1 . 329E-03 

1152 

1 . 246E-03 

1216 

1. 183E-03 

1280 

1. 119E-03 

1280 

1.913E-06 
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Plotting  for  day: 


10.00  has  been  done 


2  restart  records  written 

A  binary  restart/history  file  is  also  produced.  It  contains  the  model  fields  at  times  0  and 
10  days  as  well  as  all  the  other  information  needed  by  the  model  to  continue  the  run  where 
it  left  off. 

The  final  output  is  an  NCAR  graphics  metafile.  The  pictures  from  this  file  are  shown  in 
figures  5.2  to  5.18. 

5.3  Changes  needed  for  an  open  domain 

5.3.1  mud2 

First  of  all,  the  number  of  grid  points  is  different  for  a  periodic  domain,  as  described  in 
section  5.1.1.  The  L  dimension  must  be  one  less  to  be  consistent  with  equation  (5.1).  The 
size  of  nwrk  can  also  be  made  smaller.  For  a  periodic  channel 

nwrk  =  16  •  4  ■  L  •  M/3 

while  for  an  open  domain 

nwrk  =  14  •  4  •  L  •  M/3 
is  sufficient.  See  the  mud2  documentation. 

5.3.2  Block  data 

The  variable  perchan  should  be  set  to  false  in  the  blockdata. 

5.3.3  Boundary  conditions 

The  subroutine  bcs  must  be  changed  to  provide  the  desired  boundary  conditions.  The 
default  is  a  closed  basin  when  perchan  is  false. 

5.3.4  ezgrid 

In  ezgrid  the  line 

dx  =  xl  /  float (Lm2) 
should  be  changed  to 

dx  =  xl  /  float(Lm). 

In  ezgrid  the  variable  xl  is  the  length  of  the  periodic  domain  which  goes  from  £  =  1  to  £  = 
Lm.  in  the  SPEM  the  variable  xl  is  the  length  of  the  channel  from  £  =  1  to  <f  =  L  and  is 
only  used  to  scale  the  plots. 
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Sedmount  txdmoie 
50TTCM  ^CPOGRAP^Y 
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Figure  5.2:  Bottom  topography 
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Sedmount  Exdmole 
BOTTOM  SLOPE 
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Seamount  Examoie  MAX  RATIO  = 

H/V  RESOLUTION 


CONTOUR  FROM  .2  ”  3.-  5 


Figure  5.4:  Resolution  factor 
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Seamount  Example 

!  I  DA v  =  0.0) 


CONSTANT  FIELD  -  VALUE  IS  0 


Figure  5.5:  ip  field  at  day  0 
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Sedmount  Example 

VELCCITY (Z  =  -400.3. DAY  = 


o 


\ 


Figure  5.6:  v  field  at  day  0 
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CONTOUR  FROM  -.00006525  TO  -.00005075  BT  .C0000C5 


S(Z  =  -400.0. DAY  =  0.0)  T(Z  =  -400.0. DAY  =  0.0) 


CONTOUR  FROM  1  TO  31  BY  1  CONTOUR  FROM  3.822935  TO  3.82304  BY  . C00CC5 


Figure  5.7:  Slabs  of  the  w,p,  S  and  T  fields  at  day  0 
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Seamount  Exdmole 


StXI  =  21. DAY  =  0.0)  T (XI  =  21. DAY  =  0.0) 


CONTOUR  FROM  1  TO  31  8Y  1  CONTOUR  FROM  1 .6  TO  4.8  BY  .2 

RHO ( XI  =  21 .DAY  =  0.0) 


CONTOUR  FROM  -.000185  TO  -.000005  8Y  .00001 


Figure  5.9:  Constant  £  slices  of  the  5,  T  and  p  fields  at  day  0 


Sedmount  Example 


UiETA  -  21. DAY  =  0.0)  V1ETA  =  21. DAY  =  O.C 


UlETA  =  21 .DAY  =  0.0) 


Figure  5.10:  Constant  t]  slices  of  the  u.v  and  w  fields  at  day  0 
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Seamount  Example 


SlETA  =  21. DAY  =  0.0)  T1ETA  =  21, DAY  =  0.0) 


CONTOUR  FROM  1  .6  TO  d.8  ST  .2 


RHOtETA  =  21. DAY  =  0.0) 


CONTOUR  FROM  -.000185  TO  -.000005  8T  .00001 


Figure  5.11:  Constant  77  slices  of  the  S,  T  and  p  fields  at  day  0 
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CONTOUR  FROM  -287500  TO  237500  BY  250CC 


Figure  5.12:  ip  field  at  day  10 


Seamount  Example 

VELOCITY ; 2  -  -400,0, day  =  10.0) 


MAXIMUM  VtCTQO 


Figure  5.13:  v  field  at  day  10 
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S(Z  =  -400.0. DAY  =  10.01  T(Z  =  -400.0. DAY  =  iQ.O! 


CONTOUR  RROfl  1  TO  31  er  1  CONTOUR  RROfl  3.809  TC  3.336  3y  .CO- 


Figure  5.14:  Slabs  of  the  w,p,S  and  T  fields  at  day  10 
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Seamount  Example 


S(XI  =  21 . DAY  =  10.01  T ( X I  =  21. DAY  =  10.0! 


RHOIXI  =  21. DAY  =  10.0) 


Figure  5.16:  Constant  £  slices  of  the  S.T  and  p  fields  at  day  10 
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Seamount  Example 


IHETA  =  21  .DAY  =  10.0) 


CONTOUR  FROM  -.05125  TO  -.00125  ST  .0025 


CONTOUR  FROM  -  00085  TO  .00125  BY  .0001 


CONTOUR  TROM  -.0045  TO  .0195  BY 


Figure  5.17:  Constant  rj  slices  of  the  u,v  and  w  fields  at  day  10 


Seamount  Example 


S ( ET  A  =  21. DAY  =  10.0)  T ( ETA  =  21. DAY  =  10.0) 


CONTOUR  FROM  16.295  TO  16.375  BY  .005  CONTOUR  FROM  1 .6  TO  4.0  BY  2 

RH01ETA  =  21. DAY  =  10.0) 


CONTOUR  FROM  -.003125  TO  .001375  SY  .00025 


Figure  5.18:  Constant  77  slices  of  the  S.T  and  p  fields  at  day  10 
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Chapter  6 

Lagrangian  Floats 


6.1  General  description 

SPEM  has  optional  Lagrangian  floats  for  producing  Lagrangian  statistics  and  for  simulating 
drifters.  The  float  information  is  stored  in  an  array  and  the  floats  are  advected  at  each 
timestep  by  the  model  velocity  fields.  The  float  positions  and  values  of  T ,  S ,  p,  and  u  at 
these  positions  are  written  out  to  a  float  history  file  at  user-specified  intervals.  Plots  and 
float  statistics  can  be  created  from  this  history  file. 

The  particles  are  advected  using  a  fourth  order  Runge-Kutta  scheme.  For  this  the  fluid 
velocity  at  the  particle  position  is  needed.  The  particles  can  be  located  between  the  grid 
points  so  a  horizontal  interpolation  is  done  using  the  neighboring  grid  points,  there  is  a 
choice  of  a  bilinear  interpolation  using  the  four  nearest  grid  points  or  a  bicubic  interpolation 
using  the  four  by  four  square  of  grid  points  surrounding  the  particle.  The  algorithm  for  the 
bicubic  interpolation  is 


A{xj,y}) 


E  E 


A{xi,yj) 


nk&A*f  -  xk)  Ulit\A(yf  ~  yt) 

njbfif4(s»  -  **)  ni=Ji,4(yj  -  yi) 


(6.1) 


where  (z/,  yj)  is  the  position  of  the  particle  and  ( z;,  )  are  the  positions  of  the  16  neighboring 

grid  points.  The  interpolation  is  done  in  the  horizontal  plane  at  the  depth,  in  sigma  space, 
of  the  particle.  The  velocities  at  that  depth  are  calculated  using  the  polynomial  modes  of 
the  model. 

At  each  timestep,  each  particle  is  checked  to  see  if  it  is  inside  the  model  domain.  If  a 
particle  leaves,  it  is  reinitialized  at  a  fixed  location  which  is  specified  by  the  user  (at  xnew. 
ynew,  and  znew)  and  a  message  is  printed  to  the  standard  output  (unit  6).  If  you  find  that 
particles  are  often  leaving  the  box  other  than  through  open  boundaries,  then  your  timestep 
may  be  too  long. 


6.1.1  bcw 

In  the  model  the  vertical  velocity  array  w  ends  one  half  grid  point  to  the  inside  of  the  edge 
of  the  model  domain.  To  allow  the  same  interpolation  scheme  in  all  of  the  velocity  fields  it 
was  decided  to  add  exterior  points  to  the  w  array.  The  exterior  values  are  determined  in 
the  subroutine  bcw  by  a  linear  extrapolation  from  the  interior  values  of  the  array.  These 
exterior  values  are  used  only  in  the  particle  tracking  subroutines. 


6.2  To  use  floats 

The  changes  which  need  to  be  made  to  the  code  to  use  floats  are  as  follows: 


6.2.1  nflt 

Decide  how  many  floats  you  wish  to  follow  and  set  the  parameter  nflt  accordingly. 

6.2.2  blockdata 

In  the  blockdata,  set  the  variable  fltflag  to  true.  Make  sure  the  variables  flt4flag,  fltrrec, 
and  kptime  are  given  appropriate  values. 


6.2.3  fltinit 

Much  like  init,  fltinit  will  either  initialize  the  floats  from  a  history  file  or  it  will  provide 
initial  positions  for  ail  the  floats.  These  float  positions  are  in  the  mapped  model  coordinate 
system,  i.e.  £,77, <7.  The  floats  must  lie  in  the  range  1  <  £  <  L,  1  <  77  <  M,  and  — 1  <  a  <  1. 

Some  aspects  of  the  float  plotting  program  fltplt  assume  that  the  floats  were  released 
in  logical  groups,  for  instance,  in  five  groups  deployed  at  five  2  levels  in  the  same  pattern. 
However,  it  is  not  necessary  to  group  the  floats  in  this  way. 

6.2.4  fltstp 

In  the  subroutine  fltstp,  there  are  the  local  variables  xnew,  ynew,  and  znew.  They  specify 
the  position  at  which  a  float  will  be  reinitialized  if  it  leaves  the  model  domain.  Once  again, 
these  values  are  in  the  model  coordinate  system  (£,  77,  cr).  A  more  sophisticated  scheme  could 
be  used  in  an  open  domain  where  you  expect  the  floats  to  have  a  larger  chance  of  leaving. 

6.2.5  rhoflt 

The  subroutine  rhoflt  computes  the  density  at  the  floats  from  the  temperature  and  salinity 
at  the  floats.  Check  to  make  sure  it  is  consistent  with  rhocal. 


6.3  Plotting  the  float  tracks  with  fltplt 

At  any  given  time,  SPEM  knows  only  the  current  float  positions.  Therefore,  a  separate 
program,  called  fltplt,  was  written  to  plot  float  tracks  using  the  information  in  the  float 
history  file.  It  makes  two-dimensional  plots  in  x,  y  coordinates  and  draws  the  float  tracks  in 
color  where  color  represents  one  of  the  variables  c,u;,T,  S  or  p. 

There  is  an  option  of  plotting  a  density  slab  as  a  background  to  the  ito«..f  tracks.  This 
option  requires  access  to  the  appropriate  SPEM  history  file  and  also  the  same  equation  of 
state  which  was  used  by  SPEM. 


66 


6.3.1  Changes  to  fltplt 

Like  SPEM,  fltplt  has  parameters  and  variables  which  must  be  changed  for  each  application. 

Parameters 

nflt  Number  of  floats. 

ntm  Number  of  float  records  (or  larger),  fltplt  reads  in  all  the  float  information 
which  was  written  out  by  SPEM  and  some  arrays  have  to  be  dimensioned  large 
enough.  If  you  wrote  out  float  information  4  times/day  and  ran  for  40  days,  you 
need  ntm  to  be  at  least  160. 

nlayer  One  possible  float  release  strategy  is  to  put  out  an  array  of  floats  at  each  of 
several  depths.  All  the  floats  at  one  depth  would  then  logically  belong  together 
and  are  put  in  the  same  plot  in  the  ‘default’  plots.  The  number  of  these  float 
groupings  is  given  by  the  parameter  nlayer.  nlayer  should  be  set  to  1  if  a  different 
release  strategy  is  used. 

Main  program 

iin,  iflt,  igrid  Set  the  I/O  unit  numbers. 

gridfile  For  plots  without  a  p  slab,  the  program  only  needs  to  read  in  the  grid  file  (unit 
3)  and  the  float  file  (unit  12).  If  you  want  p  slabs  you  should  have  a  history  file 
available.  This  may  have  a  header  (with  ident,  climatology,  etc).  If  gridfile  is 
true  then  fltplt  will  try  to  read  the  grid  file,  otherwise  it  will  read  in  the  header. 

perchan,  eqstflag  Same  as  in  SPEM.  eqstflag  is  only  needed  when  plotting  p  slabs. 

fpday  Number  of  float  records  per  day.  Calculate  it  from  the  model  timestep  dt  and 
the  frequency  of  saving  float  records  kptime.  That  is. 

„  ,  S6400 

fpday  =  — - - — — - . 

dt  x  kptime 

Character  quality  The  quality  of  the  characters  diawn  by  plotchar  is  set  by  a  call 
to  pcseti.  The  best  is  0.  the  cheapest  is  2. 

lwidth  Sets  the  width  of  the  outline  drawn  by  pltperim  and  the  width  of  the  float 
tracks.  Makes  no  difference  on  some  hardware. 

rhocal  Make  sure  that  you  have  the  same  rhocal  as  in  SPEM. 

Block  data  Initially,  just  make  sure  there  are  the  proper  number  of  values  (nlayer)  for 
each  limit  array.  If  you  want  to  set  the  color  range  on  the  plots,  the  limits  can  be  set 
in  the  block  data.  To  use  vour  limits,  type  ‘iT  in  response  to  the  question  ‘Should  the 
range  of  colors  be  calculated  from  the  chosen  floats?'. 

Colors  Each  of  the  five  plots  (z,  w,T,  S, p)  has  its  own  color  bar.  defined  in  tzplt,  twplt, 
ttplt,  tsplt  and  trhoplt.  To  change  a  color  bar.  set  the  parameter  ncol  in  the 
appropriate  subroutine  to  the  number  of  colors  you  want,  then  fill  the  rgb  array. 


(Note  that  some  compilers  have  a  default  limit  to  the  number  of  allowed  continuation 
lines). 

The  foreground  and  background  colors  have  been  set  everywhere  to  black  and  white 
respectively.  Change  the  calls  to  gscr  if  you  would  like  them  to  be  different. 

6.3.2  Running  fltplt 

It  will  interactively  ask  you  a  number  of  questions,  most  of  which  should  be  obvious.  Some 
which  aren’t  are: 

•  ‘Would  you  like  the  default  plots?’  Answering  ‘yes’  will  tell  it  to  plot  all  the  floats  in 
each  group  together.  For  each  of  the  nlayer  groups,  there  will  be  a  plot  for  each  of 
z,w,T,S  and  p.  There  will  be  no  density  slabs  and  no  more  questions. 

•  ‘Would  you  like  to  specify  a  range  of  float  IDs?’  This  is  so  you  can  ask  for  floats  121 
to  140  without  typing  in  121,  122,  123,  ...  It  will  assume  that  all  floats  to  be  plotted 
have  consecutive  ID’s. 

•  ‘Should  the  range  of  colors  be  calculated  from  the  chosen  floats?’  If  not.  the  limits  set 
in  the  blockdata  will  be  used.  The  limits  will  be  that  of  the  group  to  which  the  first 
float  to  be  plotted  belongs. 

•  'Would  you  like  to  plot  a  density  surface?’  If  so,  you  should  know  the  structure  of  your 
model  field  history  file.  For  example,  if  you  want  the  density  field  from  day  60,  you 
have  to  know  that  day  60  is  in  the  seventh  restart  record. 

The  depth  is  given  as  a  2  value  which  is  a  negative  number. 

For  the  contour  interval,  you  can  say  .3  to  get  a  contour  interval  of  .3,  0.  to  accept  the 
conpack  default  of  about  16  contours,  or  a  negative  integer  ncl  to  get  about  — ncl 
contours. 

6.4  Seamount  example 

To  demonstrate  the  use  of  floats,  the  process  of  putting  them  in  the  seamount  resonance 
version  from  section  5.2  will  be  described.  However,  this  application  will  not  lead  to  the 
most  dramatic  float  tracks  since  the  flow  speeds  are  so  small. 

6.4.1  In  SPEM 

nflt  Let’s  put  a  5  by  5  array  of  floats  over  the  seamount  at  two  depths,  using  50  floats. 

blockdata  Set  fltflag  to  true  and  flt4flag  to  false  (for  economy),  fltrrec  is  0  initially. 
To  produce  8  float  records  per  day.  set  kptime  to  16. 

fltinit  The  center  of  the  seamount  is  if  (£,  //)  =  (21.5,  21.5),  which  will  also  be  the  center 
of  the  float  array.  The  float  spacing  will  be  A£  =  A/?  =  2.  The  first  group  of  25  floats 
will  be  at  a  depth  of  —100  meters  while  the  second  group  will  be  at  —499.9  meters. 


fltstp  The  floats  should  not  be  leaving  the  domain  so  the  restart  location  (xnew.ynew.znew) 
is  not  critically  important.  Floats  will  get  restarted  in  the  corner  at  (£,  77,  a)  =  (2, 2,  .9). 


rhoflt  The  equation  of  state  in  rhoflt  is  consistent  with  that  in  rhocal. 

6.4.2  In  fltplt 

Parameters  Set  the  parameters  nflt,  ntm  and  nlayer  to  50,  81  and  2  respectively. 

Main  program  We  want  to  plot  density  slices  with  the  float  tracks  so  we  will  be  reading 
the  SPEM  history  file.  It  has  a  header  with  the  grid  information  so  gridfile  is  false. 
The  SPEM  history  file  unit  is  specified  by  iin  and  is  set  to  1.  The  float  history  file 
unit  is  specified  by  iflt  and  is  set  to  12. 

perchan  and  eqstflag  are  true  and  false  respectively,  fpday  is  8. 

rhocal  Use  the  same  equation  of  state  as  in  SPEM  for  the  density  slab  plots. 

Block  data  Make  sure  there  are  two  values  for  each  of  the  limits  arrays,  one  for  each  float 
group. 

Colors  The  color  bars  each  have  14  colors,  7  reds  and  7  blues. 

6.4.3  Example  plots 

There  are  many  options  when  running  fltplt.  To  get  an  overview  of  what  all  the  floats  do.  it 
is  convenient  to  ask  for  the  default  plots,  the  first  frame  of  which  is  shown  in  figure  6.1.  It 
is  then  possible  to  focus  in  on  a  few  floats  and  to  overlay  a  density  slab  plot  to  get  an  idea 
of  how  the  float  motion  relates  to  the  physical  fields,  as  shown  in  figure  6.2. 
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Mini  mum  Z  -  -101.2 

Maximum  Z  -  -99.4 


Min  Delta 
Max  Delta 


Seamount  Examoie 
Float  IracKs  for  Days 
0.0  to  10.0 


Figure  6.1:  First  frame  of  the  default  plots 
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Density  Field  from  Ddy  10.0  Depth  -pQO.Q 


CONTOUR  FROM  27  694498  TO  27.698997  BY  .000499 


Figure  6.2:  Several  floats  and  a  density  surface 


Appendix  A 
Model  Time-step 


Numerical  time-stepping  uses  a  discrete  approximation  to 

du(t) 


dt 


=  m 


(A.l) 


where  F(t)  represents  all  the  right-hand  side  terms.  The  simplest  approximation  is  the  Euler 
time  step: 

u(t  +  A£)  —  u(t) 


A  t 


=  T{t) 


(A. 2) 


where  you  predict  the  next  u  value  based  only  on  the  current  fields.  This  method  is  accurate 
to  first  order  in  A t  and  is  stable  for  a  sufficiently  short  time-step,  dictated  by  a  combination 
of  the  advective  speeds,  friction  coefficients,  and  the  horizontal  and  vertical  grid  spacings. 
The  leapfrog  time-step  is  accurate  to  0(A t2): 


u(t.  +  At)  —  u(t  —  At) 

2A t 


(A.  3) 


This  time-step  is  more  accurate,  but  the  even  and  odd  time  steps  tend  to  diverge  in  a 
computational  mode.  This  computational  mode  can  be  damped  by  taking  an  occasional 
correction  step.  The  SPEM  has  a  choice  of  correction  steps,  both  of  which  use  a  leapfrog 
step  to  obtain  an  initial  guess  of  u{t  +  At).  We  will  call  the  right-hand  side  terms  calculated 
from  this  initial  guess  tT'(t  +  At). 

The  first  choice  of  correction  steps  is  a  trapezoidal  step: 


«( t  +  At)  —  u(t)  1  . 

- - Yt - -  =  ^  [F(t)  +  f'(t  +  At)] 

and  the  second  choice  is  a  third-order  Adams- Bashforth: 
u( /  -f-  At )  —  //( t ) 


At 


=  T{t)  -  -7(t  -  At)  +  -T’{t  +  A r 


(A.l) 


A. 5 ) 


The  SPEM  uses  a  leapfrog  time  step  with  a  correction  every  ncorr  time  steps.  The  type 
of  correction  step  is  determined  by  the  value  of  stpflag  (true  for  a  third-order  Adams- 
13  ashforth). 

The  model  carries  two  time  levels  for  many  of  the  physical  fields,  differentiated  by  the 
last  index  in  the  arrays.  The  initial  fields  are  mily  read  in  for  one  time  level:  an  Euler  time 
step  is  used  for  the  first  step  to  get  things  going.  A  correction  step  is  made  if  you  are  using 
trapezoidal  correct  ions  (you  would  need  to  obtain  J-{  —  A/ )  to  make  a  t  bird  order  correct  ion  h 

Every  nrst  time  steps  a  restart  record  is  written.  If  you  were  to  start  the  model  from  one 
of  these  restart  records  you  would  start  with  an  Euler  step.  Therefore,  to  give  repeatable 
result >.  the  model  treats  the  first  -tep  after  each  record  ;s  written  as  it  it  were  the  first  step. 


Appendix  B 

Chebyshev  Polynomial  Basis 
Functions 


As  described  in  section  3.1,  the  P*(<x)  expansion  functions  used  in  the  vertical  spectral 
method  can  be  chosen  somewhat  arbitrarily.  Here  we  review  the  form  of  the  matrix  ooerators 
when  the  Pfc(cr)  are  a  modified  form  of  the  Chebyshev  polynomials  of  the  first  kind  (Tt(rr): 
see,  e.g.,  [1]  and  [9]).  The  T^(o)  are  defined  as 

Tfc(cr)  =  cos|P  cos-1  (cr)],  —  1  <  o  <  1 

where  —it  <  cos~l(a)  <  0.  We  then  set 


Pfc(cr)  = 


Ua) 

o 

II 

TM 

k  >  1,  k  odd 

(B.l) 

Tk{<?)  +  fc2_1 

k  >  2,  k  even. 

which  then  conform  to  the  required  integral  property  that  only  P0(a)  have  a  non-zero  vertical 
integral.  The  collocation  points  an  are  located  at  the  extrema  of  the  highest  order  polynomial 
P  v  ( <7 ) ;  hence. 

crn  =  cos[7r(n  —  N)/N],  0  <  n  <  N.  ( B . 2 ) 

The  first  eight  modified  Chebyshev  polynomials  and  the  collocation  levels  for  .V  =  8  are 
shown  in  figure  B.l. 

The  matrix  F  (see  equation  (3.3))  may  now  be  evaluated  form  (B.l)  and  (B.2).  Substi¬ 
tuting  (B.l)  into  (3.6)  gives  the  elements  of  matrix  R: 


R-nk  = 


OPk(cr)  ks\n(kOn) 


do- 


sin  0n 


(B.3) 


where  0n  —  cos  x{on).  Similarly.  (B.l)  and  (3.8)  give  the  elements  of  matrix  5: 


Sak  = 


Pt((T  )do  =  < 


(1  ~<rn) 

§(1-0 

-os(  k  +  1  )'7  cost  1  )rr  1  ^ 

2(^  +  1)  2(k- 1)  J,;os-l(<,n) 

josffc+ll'T  ros(t— 11^  ^ 

lk+Ti  2[k-\)  J, 

+(1  -  cr,P  1 


k‘-l 


k  =  0 
k  =  1 

k  >  l.odd 

k  >  1 . even. 


(B.l) 


f  inally,  the  differentiating  and  integrating  operators  Co/  and  (’isr  are  computed  by  first 
inverting  F  and  calculating  the  products  (  ’ o/  =  f)F~l  and  C[ \j  =  FF~l. 


Appendix  C 


Topographic  Steepness 


All  sigma  coordinate  models  are  known  to  have  some  trouble  in  the  presence  of  steep  to¬ 
pography,  especially  if  the  grid  spacing  is  not  fine  enough.  Unfortunately,  we  do  not  know 
in  general  the  relationship  between  steep  topography,  stratification,  horizontal  and  vertical 
resolution  and  whether  the  model  will  run  stably.  However,  two  plots  have  been  added 
to  ploth  to  guide  the  user  in  determining  which  areas  of  the  domain  are  likely  to  lead  to 
topographic  problems. 

The  first  of  these  plots  is  simply  the  bottom  slope,  or  |V/j|.  The  maximum  for  the  whole 
domain  is  written  in  the  corner  of  the  plot.  We  have  successfully  run  the  SPEM  with  bottom 
slopes  of  up  to  ~  0.10,  while  continental  shelves  can  have  slopes  in  excess  of  0.25. 

If  you  find  that  the  model  is  ill-behaved  and  you  suspect  a  problem  with  steep  topography, 
there  are  two  ways  to  improve  on  the  situation.  One  is  to  filter  the  topography  until  it  is  no 
longer  ‘too  steep’.  The  other  is  to  increase  the  horizontal  grid  spacing  until  you  resolve  the 
slope  ‘well  enough’.  A  measure  of  how  well  the  slope  has  to  be  resolved  is  gotten  from  the 
“7£-value”  (the  ratio  of  horizontal  to  vertical  resolution): 

Ah 

A(hcr) 

This  is  a  three-dimensional  quantity,  but  we  only  plot  it  as  a  function  of  x  and  y  since 
we  know  that  when  using  the  Chebyshev  polynomials,  the  maximum  is  between  the  two 
collocation  levels  closest  to  the  bottom. 

7 Z  is  related  to  the  "hydrostatic  consistency”  requirement  known  from  discrete  level 
models,  where  TZ  has  to  be  less  than  1,  but  SPEM  has  run  with  values  up  to  5.  Once  again, 
we  cannot  tell  exactly  what  will  or  will  not  run.  These  tools  are  provided  to  help  diagnose 
problems  with  steep  topography. 


Appendix  D 

Context  Diffs  for  velvet 


What  follows  is  a  Unix  context  difference  file.  It  starts  with  a  header  specifying  which  files  are 
being  compared.  The  lines  which  differ  are  shown  with  a  few  surrounding  lines  for  context 
and  are  marked  with  an  exclamation  point.  The  lines  from  the  distributed  velvet. f  are 
shown  first,  followed  by  the  lines  from  the  modified  version.  See  section  4.5  for  a  description 
of  the  changes. 

diff  -c  velvct/velvct .f  velvet .mod/velvct . f 
***  velvct/velvct .f  Fri  Jan  26  09:24:17  1990 

-  velvet .mod/velvct .f  Mon  Mar  19  09:52:58  1990 

*************** 

***  406,416  **** 


C 

C 


ILAB  IS  NON-ZERO. 


SAVE 

FX(XX.YY)  =  XX 
FY (XX , YY)  =  YY 

DIST(XX.YY)  =  SQRT(XX*XX+YY*YY) 

MXF (XX , YY , UU , W , SFXX , SFYY , MXX , MYY )  =  MXX+IFIX(SFXX*UU) 

MYF (XX , YY , UU , VV , SFXX , SFYY , MXX , MYY)  =  MYY+IFIX(SFYY*VV) 

SCALEX (MM , NN , INCXX , INCYY , HAA , XX 1 , XX2 , YY 1 , YY2 , XX3 , XX4 , YY3 , YY4 , 
1  LENN)  =  LENN/HAA 

SCALEY (MM , NN , INCXX , INCYY , HAA , XX 1 , XX2 , YY 1 , YY2 , XX3 , XX4 , YY3 , YY4 , 
—  406,416  - 


C 

c 


ILAB  IS  NON-ZERO. 


c 

c 


c 

c 


SAVE 

FX(XX, YY)  =  XX 
FY (XX, YY)  =  YY 

DIST(XX, YY)  =  SQRT (XX*XX+YY*YY) 

MXF(XX,YY,UU,VV, SFXX, SFYY, MXX, MYY)  =  MXX+IFIX(SFXX*UU) 

MYF (XX, YY.UU.VV, SFXX, SFYY, MXX, MYY)  =  MYY+IFIX(SFYY*VV) 

SCALEX (MM , NN , INCXX , INCYY , HAA , XX 1 , XX2 , YY1 , YY2 , XX3 , XX4 , YY3 , YY4 , 
1  LENN)  =  LENN/HAA 

SCALEY (MM , NN , INCXX , INCYY , HAA , XX 1 , XX2 , YY 1 , YY2 , XX3 , XX4 , YY3 , YY4 , 
*************** 

***  551,556  **** 

—  551,561  - 

ZMN  =  LEN* (GL/HA) 


ZMX  =  FLOAT (LEN)+. 01 


C 

+  c  turn  off  clipping  so  vectors  can  go  out  of  domain 
+  c 

+  CALL  GQCLIP(IER,ICLP,IAR) 

+  CALL  GSCLIP(O) 

+  c 

C  DRAW  THE  VECTORS. 

C 

DO  123  J-1,NY,INCY 

***  580,587  **** 

C 

C  TURN  OFF  CLIPPING  SO  ARROW  CAN  BE  DRAWN 
C 

!  CALL  GQCLIP(IER,ICLP ,IAR) 

!  CALL  GSCLIP(O) 

CALL  DRWVEC  (28768, 608, 28768+LEN, 608, LABEL, 10) 
C 

C  RESTORE  CLIPPING 

—  585,592  - 

C 

C  TURN  OFF  CLIPPING  SO  ARROW  CAN  BE  DRAWN 

C 

!  c  CALL  GQCLIP (IER , ICLP , IAR) 

!  c  CALL  GSCLIP(O) 

CALL  DRWVEC  (28768 , 608 , 28768+LEN, 608, LABEL, 10) 
C 

C  RESTORE  CLIPPING 


SO 


L 
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